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ABSTRACT 


Two  algorithms  are  developed  which  address  the  problem  of  estimat- 
ing the  magnitude  and  phase  of  the  optical  transfer  function  associated 
with  a blurred  image.  The  primary  focus  of  the  research  is  on  the 
estimate  of  the  phase  of  the  optical  transfer  function.  With  the 
sharpening  of  one  approximation,  the  method  affords  a reasonable 
estimate  of  the  phase  of  the  optical  transfer  function.  Once  an 
estimate  of  the  optical  transfer  function  has  been  made,  the  correspond- 
ing blurred  image  is  Wiener  filtered  to  estimate  the  original  unblurred 
image.  Results  are  demonstrated  on  computer  simulated  blurs  and  also 
on  real  world  blurred  imagery.  Included  is  a mathematical  bound  on 
the  phase  of  the  optical  transfer  function. 


I 


i ' 


V 


j 


TABLE  OF  CONTENTS 


Chapter 


INTRODUCTION 

1.1  Introduction  and  Literature  Review 

1.2  Research  Goals 

THE  ALGORITHMS  TO  BE  INVESTIGATED 
NATURE  OF  THE  PHASE 

3.1  Discussion 

3.2  Phase  Bound 
RESTORATION 

IMPLEMENTATION  AND  SYSTEM  PERFORMANCE 

5.1  Introduction 

5.2  On  the  Assumed  Model 

5.3  On  the  Question  of  Estimating  Statistical 
Properties  of  the  Unknowns/Convergence 

5.4  On  the  Stability  of  the  Iterations 

5.5  On  Errors  of  Approximation 

5.6  On  Phase  Averaging 

5.7  On  the  Relaxation  of  One  Assumption 

5.8  Summary 

RESULTS  OF  SIMULATIONS 

6.1  Presentation  of  Results 

6.2  Discussion 

RESULTS  ON  REAL  WORLD  BLURRED  IMAGES 


SUMiARY 


vi 


Chapter 


TABLE  OF  CONTENTS  (CONT'D) 


Page 


8.1  Suirmary  187 

8.2  Open  Question  I87 

Appendix  A.  On  a Theorem  Concerning  the  Zeros  of  a 

Trigonometric  Polynomial  189 

REFERENCES  192 


‘ 1 


I 


I 


r 


LIST  OF  FIGURES 


Figure 

Page 

3.1 

Principal  value  of  the  phase  of  the  Fourier 
transform  of  a 64  x 64  pixel  subimage 

17 

3.2 

Phases  corresponding  to  two  functions 

18 

3.3 

Phases  corresponding  to  two  functions 

19 

3.4 

Phases  corresponding  to  two  functions 

20 

3.5 

Phases  corresponding  to  two  functions 

21 

3.6 

Function  of  Figure  3.2  together  with  functions 
corresponding  to  10%,  20%,  30%,  40%,  and  50% 
phase  distortions 

23 

3.7 

Function  of  Figure  3.3  together  with  functions 
corresponding  to  10%,  20%,  30%,  40%,  and  50% 
phase  distortions 

24 

3.8 

Function  of  Figure  3.4  together  with  functions 
corresponding  to  10%,  20%,  30%,  40%,  and  50% 
phase  distortions 

25 

3.9 

Function  of  Figure  3.5  together  with  functions 
corresponding  to  10%,  20%,  30%,  40%,  and  50% 
phase  distortions 

26 

3.10 

Original  image 

28 

3.11 

Phase  distorted  image 

28 

3.12 

Phase  distortion 

28 

3.13 

A once  sharp  image  reduced  to  garbage 

29 

3.14 

Location  of  sample  points 

31 

3.15 

Representation  of  H(u) 

32 

3.16 

Representation  of  H(u)  with  R(u)  zero  on  an  interval 

32 

3.17 

Example  of  phase  bound  for  L=3 

36 

4.1 

Schematic  of  restoration  process  using  interpolation 
in  the  Fourier  domain 

41 

4.2 

Restoration  using  64  x 64  pixel  sub-blocks 

42  Vi 

tAi 


LIST  OF  FIGURES  (CONI' D) 


Figure 

Page 

4.3 

Restoration  using  approach  diagrammed  in 

Figure  4.1 

43 

4.4 

Triangularly  shaped  PSF 

45 

4.5 

Comparison  of  approaches  to  restoration 

46 

4.6 

Comparison  of  approaches  to  restoration 

47 

5.1 

Four  statistically  similar  images 

51 

5.2 

Magnitude  autocorrelation  for  Au=l,  Av=0 

53 

5.3 

Magnitude  autocorrelation  for  Au=0,  Av=l 

54 

5.4 

Phase  difference  histogram  for  case  of  no 
windowing,  au=1,  av=0 

57-58 

5.5 

Phase  difference  histogram  for  case  of  no 
windowing,  au=0,  av=1 

59-60 

5.6 

Phase  difference  histograms  for  case  of  no 
windowing,  on  axis,  Au=l ,Av=0 

63-64 

5.7 

Phase  difference  histograms  for  case  of  no 
windowing, on  axis,  Au=0,  Av=l 

65-66 

5.8 

Phase  difference  histograms  for  case  of  Parzen 
windowing,  au=1 , Av=0 

67-68 

5.9 

Phase  difference  histograms  for  case  of  Parzen 
windowing,  Au=0,  Av=l 

69-70 

5.10 

Standard  deviation  of  the  average  vs.  number  of 
samples  for  a sample  standard  deviation  of  54° 

72 

5.11 

Comparison  of  magnitude  estimates  using  knowledge  of 
the  undegraded  image 

77 

5.12 

Restorations 

78 

5.13 

Comparison  of  magnitude  estimates  without  knowledge 
of  the  undegraded  image 

79 

5.14 

Phases  corresponding  to  a function  convolved 
aperiodically  and  circularly 

86 

LIST  OF  FIGURES  (CONT'D) 


Page 


Figure 


Phases  corresponding  to  a function  convolved 
aperiodically  and  circularly  87 

Phases  corresponding  to  a function  convolved 
aperiodically  and  circularly  88 

Phases  corresponding  to  a function  convolved 
aperiodically  and  circularly  89 

Candidate  windows  90 

Arbitrarily  chosen  subimage  used  in  windowing  study  91 

Magnitude  differences  of  Fourier  transforms  of 
subimage  convolved  aperiodically  and  circularly  for 
different  windows  92 

Phase  differences  of  Fourier  transforms  of  subimage 
convolved  aperiodically  and  circularly  for  different 
windows  94 

Simulation  results  corresponding  to  no  windowing  99 

Simulation  results  corresponding  to  no  windowing  100 

Simulation  results  corresponding  to  Parzen  windowing  102 

Simulation  results  corresponding  to  Parzen  windowing  103 

Vector  sum  A + B 108 

Average  vector  B corresponding  to  H(I)  108 

The  resultants  A^+B^  109 

Representation  of  the  complex  value  H(u^,v^)  111 

Vectors  of  average  phase  of  roughly  it  radians  111 

Phase  differences  corresponding  to  OTF,  undegraded  114 

subimages,  and  degraded  subimages 


5.33 


Phase  zones 

An  image  zeroed  at  the  edges  of  the  subimages  and 
then  bl urred 


X 


LIST  OF  FIGURES  (CONI' D) 


Figure 

Page 

5.34 

Corresponding  subimage  to  which  average  phase 
differences  of  zeroed  Image  converge  to 

119 

5.35 

Phase  plane 

121-122 

6.1 

PSF  corresponding  to  motion  blur 

131 

6.2 

Comparison  of  magnitude  of  OTF  and  estimate 

132 

6.3 

Comparison  of  phase  of  OTF  and  estl'nates 

133 

6.4 

Comparison  of  phase  of  OTF  and  estimates 

134 

6.5 

Restorations 

135-136 

6.6 

PSF  of  square  blur 

137 

6.7 

Comparison  of  magnitude  of  OTF  and  estimate 

138 

6.8 

Comparison  of  phase  of  OTF  and  estimates 

139 

6.9 

Comparison  of  phase  of  OTF  and  estimates 

140 

6.10 

Restorations 

141-142 

6.11 

PSF  corresponding  to  double  exposure 

143 

6.12 

Comparison  of  magnitude  of  OTF  and  estimate 

144 

6.13 

Comparison  of  phase  of  OTF  and  estimates 

145 

6.14 

Comparison  of  phase  of  OTF  and  estimates 

146 

6.15 

Restorations 

147-148 

6.16 

PSF  corresponding  to  quadruple  exposure 

14*9 

6.17 

Comparison  of  magnitude  of  OTF  and  estimate 

1 50 

6.18 

Comparison  of  phase  of  OTF  and  estimates 

151 

6.19 

Comparison  of  phase  of  OTF  and  estimates 

152 

6.20 

Restorations 

153-154 

6.21 

Triangularly  shaped  PSF 

155 

x1 

k — — ■ - 


LIST  OF  FIGURES  (CONI' D) 


Figure 

Page 

6.22 

Comparison  of  magnitude  of  OTF  and  estimate 

156 

6.23 

Comparison  of  phase  of  OTF  and  estimates 

157 

6.24 

Comparison  of  phase  of  OTF  and  estimates 

158 

6.25 

Restorations 

159-160 

6.26 

Triangularly  shaped  PSF 

161 

6.27 

Comparison  of  magnitude  of  OTF  and  estimate 

162 

6.28 

Comparison  of  phase  of  OTF  and  estimates 

163 

6.29 

Comparison  of  phase  of  OTF  and  estimates 

164 

6.30 

Restorations 

165-166 

7.1 

Scene  before  and  after  photographically  induced  blur 

168 

7.2 

Estimates  of  magnitude  and  phase  of  OTF 

169 

7.3 

Blurred  image  and  restorations 

170 

7.4 

Blurred  image  and  restorations 

171 

7.5 

Blurred  image  and  restorations 

172 

7.6 

Blurred  image  and  restorations 

173 

7.7 

Scene  before  and  after  photographically  induced  blur 

175 

7.8 

Estimates  of  magnitude  and  phase  of  OTF 

176 

7.9 

Blurred  image  and  restorations 

177 

7.10 

Blurred  image  and  restorations 

178 

7.11 

Blurred  image  and  restorations 

179 

7.12 

Blurred  image  and  restorations 

180 

7.13 

Estimate  of  magnitude  and  phase  of  OTF 

182 

7.14 

Blurred  image  and  restorations 

183 

7.15 

Blurred  image  and  restorations 

184 

xii 


Chapter  1 
INTRODUCTION 


1.1  Introduction  and  Literature  Review 

The  emergence  of  computer  technology  has  drastically  changed 
methodologies  In  diverse  fields.  As  the  cost  per  calculation  has  been 
reduced  by  orders  of  magnitude,  techniques  which  were  once  considered 
too  costly  are  now  economically  feasible.  In  addition,  concomitant 
technologies  have  emerged.  One  of  these  concomitant  technologies  has 
come  to  be  known  as  "digital  Image  processing."  This  dissertation  Is 
concerned  with  an  area  of  digital  Image  processing  termed  Image 
restoration;  simply  put,  the  removal  of  a blur  or  degradation  from  an 
Image.  Everyday  examples  of  blurred  Imagery  would  Include  photographs 
which  were  taken  with  an  out  of  focus  camera  or  motion  blur,  photographs 
which  were  taken  while  the  camera  and/or  the  object  were  moving. 

Although  blurred  photographs  are  the  obvious  examples,  the 
technology  of  Image  restoration  can  be  applied  to  any  two  dimensional 
display  of  data,  for  example,  radar  maps,  sonar  maps,  and  chest  x-rays. 

A mathematical  model  of  the  degrading  process  Is  as  follows.  Let 
f(x,y}  denote  the  undegraded  Image,  g(x,y}  denote  the  degraded  or 
blurred  Image,  and  n(x,y)  denote  additive  noise.  The  blurring  process 
may  be  modeled  by  equation  (1-1). 

00 

g(x,y)  =jjj|'  h(x,y,c,n)f(e.n)dcdn  + n(x,y)  (1-1) 

• CD 

The  function  h(x,y,c,n)  Is  termed  the  point  spread  function  (PSF).  Note 
that  given  the  degraded  Image  g(x,y),  one  desires  to  obtain  the 
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undegraded  Image  f(x,y}. 
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A simpler  form  of  the  PSF,  h(x,y,c,n),  is  where  h is  a function 
only  of  the  differences  between  respective  coordinates.  That  is, 

h(x,y,C,n)  = h(x-C,y-n)  . 

When  this  is  the  case,  the  PSF  is  said  to  be  spatially  invariant.  The 
significance  of  the  PSF  being  spatially  invariant  lies  in  the  fact  that 
the  blur  is  unchanged  across  the  image.  In  addition,  the  degraded  and 
undegraded  images  are  related  via  convolution: 


9(x,y)  h(x-C,y-n)f(C,n)dcdn  + n(x,y)  . 

Or  equivalently, 

00 

g(x,y)  *j|j^  h(c,n)f(x-C,y“ri)dcdT)  + n(x,y)  . 

*00 

The  convolution  integral  may  be  Fourier  transformed  to  yield 


(1-2) 
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G(u,v)  = H(u,v)F(u,v)  + N(u,v)  (1-3) 

where  the  capital  letters  denote  the  Fourier  transforms  of  their 
respective  functions  represented  by  lower  case  letters.  Thus,  when  the 
PSF  is  spatially  invariant,  one  has  the  advantage  of  being  able  to 
utilize  the  relationship  in  the  spatial  domain  (equation  (1-2))  or  the 
relationship  in  the  frequency  domain  (equation  (1-3)).  H(u,v)  will  be 
referred  to  by  the  term  optical  transfer  function  (OTF).  Note  that 
h(x,y)  and  H(u,v)  form  a Fourier  pair.  As  a result,  knowledge  of 
either  one  Implies  knowledge  of  the  other. 


Assuming  knowledge  of  h(x,y,C,n)>  one  must  solve  the  integral 
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equation  (1-1)  to  obtain  the  undegraded  Image  f(x,y).  Many  examples  of 
successful  restorations  given  a priori  knowledge  of  h(x,y.Ctn)  can  be 
found  In  the  recent  literature  [1-18]. 

In  contrast,  In  the  majority  of  practical  situations  one  would  not 
have  knowledge  of  h(x,y,c,n)>  Thus,  If  one  Is  to  utilize  techniques 
which  are  known  to  be  successful  for  estimating  f(x,y)  given  g(x,y) 
and  h(x,y,c,n),  one  must  first  estimate  h(x,y,c,n). 

To  date  techniques  for  estimating  h(x,y,c,n)  from  the  degraded 
Image  have  concentrated  on  spatially  Invariant  point  spread  functions, 
and  hereafter  h(x,y,5,n)  will  be  assumed  to  be  of  the  form  h(x-c,y-n). 
When  It  Is  clear  from  Its  context,  h(x-£,y-ii)  will  sometimes  be 
denoted  by  h(x,y). 

One  possibility  of  estimating  h(x,y)  Is  that  of  obtaining  the 
estimate  from  an  Isolated  unresolved  point  which  Is  known  to  exist 
In  the  undegraded  Image.  Accordingly,  the  Isolated  region  of  the 
undegraded  Image  can  be  modeled  as 

T(x,y)  % 6(x-t,y-n) 

where  6 denotes  a Dirac  delta  function  [19].  Thus, 

Qlsoiated^^’^^  % Jj  h(c,n)6 (x-e,y-Ti)dcdn 

••00 

^ h(x,y)  . 

This  technique  Is  obviously  quite  limited  In  that  It  assumes  the 
existence  of  an  Isolated  unresolved  point  which  Is  known  to  exist  In 
the  undegraded  Image. 

Another  possibility  Is  that  of  estimating  h(x,y)  from  sharp  edges 
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which  are  assumed  to  exist  In  the  undegraded  Image  [20].  If  h Is  one 
dimensional  and  the  edge  Is  perpendicular  to  this  dimension.  It  can  be 
shown  that 


where  g(x)  Is  defined  over  the  degraded  edge  and  the  dimension  x Is 
perpendicular  to  the  edge. 

If  one  assumes  the  possible  point  spread  functions  are  from  a 
limited  set,  one  may  use  this  knowledge  to  Idenfity  h(x,y).  For 
example,  the  spectrum  of  an  out  of  focus  image  will  have  zeros  char- 
acteristic of  an  out  of  focus  OTF.  Namely,  of 

J, (2nar) 

2 2 2 

where  denotes  a first  order  Bessel  function,  r = u +v  , and  a Is  a 
constant  related  to  the  severity  of  the  blur.  Note  that  the  zeros  will 
be  circularly  symmetric  and  functionally  related  to  parameter  a. 
Similarly,  the  spectrum  of  a motion  blurred  image  will  have  zeros 
characteristic  of  a motion  blurred  OTF.  Taking  the  spectrum  of  a 
blurred  Image  and  displaying  the  results  may  be  utilized  to  Identify 
some  blurs  [21-23]. 

An  automated  version  of  this  idea  was  demonstrated  by  Cannon 
[26-28].  He  assumed  a set  of  possible  blurs  of  out  of  focus,  motion, 
and  Gaussian,  and  automatically  via  the  computer  program  he  had  coded, 
identified  the  blurs  and  parameters  associated  with  the  blurs. 

A technique  developed  by  Stockham,  Cole,  and  Cannon  [25-28] 
which  has  been  moderately  successful  Is  as  follows.  The  degraded  image 


Is  divided  Into  subimages  which  may  overlap.  If  one  assumes  that  the 
extent  of  the  PSF  Is  small  compared  to  the  extent  of  the  subimage  and 
Ignoring  the  noise  term,  then  approximately 


g^(x,y)  !\!  I f h(x-t,y-n)f^(t,n)dcdn 

V 


(1-4) 


where  the  Index  1 denotes  the  1-th  subimage.  Taking  the  Fourier  trans- 
form of  relationship  (1-4)  one  obtains 


S^(u.v)  % H(u,v)F^(u,v)  . 

Note  that  In  the  Fourier  domain  the  functions  are  now  complex. 
Expressing  relationship  (1-5)  In  magnitude  phase  form 


(1-5) 


. . J86<(“*'')  . Jen(u.v) 

|G^(u,v)|e  ^ % |H(u,v)|e  ” 1 


jep.(u,v) 

F^(u.v)|e.  ^ (1-6) 


Cole  [25]  obtained  a reasonable  estimate  of  |H(u,v)|  as  follows.  From 
equation  (1-6) 

|G^(u,v)|  % |H(u.v)|lF^(u,v)| 

ln|6^(u,v)|  % ln|H(u,v)l  + lnlF^(u,v)|  . 

Averaging  over  the  N subimages: 

N N 

ln|G^(u,v)l  % ln|H(u,v)|  + J ln|F^(u,v)|  . 

1*1  1»1 

Assuming  the  degraded  Image  was  not  so  degraded  that  one  could  not  tell 
the  general  class  to  which  the  undegraded  Image  belonged,  It  was  shown 
experimentally  that  one  could  use  an  undegraded  Image  of  the  same  pro- 
totype class  as  the  degraded  Image  as  an  estimate  of 
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1 N 

g V 1nlF^(u.v)| 


That  Is, 


1-1  1-1 

where  denotes  subimages  of  a prototype.  Thus, 

N N 

1n|H(u,v)|  % R 2^  ■ R V 

1-1  frf 


v)|  . 


Cole's  restorations  assumed  optical  transfer  functions  of  zero  phase. 
That  Is,  e^(u,v)  - 0. 

Cannon  [26-28]  assumed  statlonarlty  and  used  the  relationship 

2 

cpg(u,v)  - ®^(u,v)|H(u.v)|  +CP^(u,v) 

to  calculate  the  magnitude  of  the  OTF  where  cp  denotes  power  spectral 
density,  cp^  was  estimated  from  g(x,y)  and  cp^  was  estimated  from  an 
undegraded  prototype  Image.  By  limiting  the  possible  blurs  to  out  of 
focus,  motion,  and  Gaussian,  Cannon  was  able  to  Identify  the  blur  and 
the  parameters  associated  with  the  given  blur.  Thus,  Cannon  was  able 
to  estimate  both  the  magnitude  and  phase  of  the  OTF  for  the  three 
possible  blurs. 

The  estimate  of  the  phase  of  the  OTF,  e^(u,v), for  the  general 
spatially  Invariant  case  has  proven  to  be  a more  difficult  problem. 
Equation  (1-6)  together  with  the  relative  success  achieved  In  estimating 
the  magnitude  of  the  OTF  suggests  a similar  approach  with  respect  to  the 
phase  of  the  OTF.  That  Is,  from  equation  (1-6) 


eg^(u,v)  % ej^(u,v)  + ep^(u,v)  . 

Averaging,  where  the  bar  denotes  averaging  in  some  sense  yet  to  be 
defined,  one  obtains 

eQ(u,v)  % e^(u,v)  + 6p(u,v)  . 

This  approach  was  considered  and/or  tried  by  a number  of  researchers 
[24,25,29].  The  basic  problem  is  that  ep(u,v)  does  not  converge  to 
anything  meaningful,  and  accordingly,  ep(u,v)  cannot  be  estimated  by  a 
prototype  or  otherwise. 

A method  which  is  closely  related  to  the  estimate  of  ej^(u,v)  is  a 
technique  studied  by  Knox  [30-32].  Knox  was  concerned  with  obtaining 
clear  photographs  of  astronomical  objects.  Since  for  earth  bound 
telescopes  astronomical  objects  must  be  viewed  through  the  atmosphere, 
the  clarity  of  these  photographs  is  generally  limited  by  atmospheric 
turbulence.  Knox  proposed  a technique  whereby  many  short  exposure 
photographs  of  the  same  object  together  with  short  exposure  photographs 
of  a nearby  star  were  taken.  These  photographs  would  then  be  digitally 
processed  to  remove  the  effects  of  the  turbulence.  Denoting  the  i-th 
short  exposure  photograph  by  g^(x,y)  and  its  corresponding  short  expos- 
ure point  spread  function  by  h^(x,y)  we  have 


g^(x,y)  « h^{x,y)*f{x,y) 


Now  consider  the  following  autocorrelation 

6^(u,v)G*(u+&u,v+Av)  = H^(u,v)H*(u+Au,v+av)F(u,v)F*(u+Au,v+av)  , 

where  the  superscript  * denotes  complex  conjugation.  Averaging  over 
the  multiple  images,  we  have 

G^(u,v)G*(u+au,v+av)  = H. (u,v)H*(u+Au,v+Av)F(u,v)F*(u+au,v+av) 

(1-8) 

where  the  bar  denotes  averaging.  Taking  the  phase  of  equation  (1-8), 
we  have 

Phase{G^ (u,v)G*(u+Au,v+Av)}  = Phase{H^(u,v)H*(u+Au,v+Av)l 

+ Phase{F(u,v)F*(u+Au,v+Av)}  (1-9) 

The  Phase  on  the  left  hand  side  of  equation  (1-9)  can  be  calculated 
from  the  multiple  images  and  the  Phase{H^(u,v)H*(u+Au,v+Av)}  can  be 
:hown  to  be  negligible  for  atmospheric  turbulence.  Denoting  the  phase 
of  F(u,v)  by  0p(u,v),  thus,  we  have 

0p(u,v)-0p(u+Au,v+Av)  % Phase{G^(u,v)G*(u+Au,v+Av)}  (1-10) 

Note  that  0p(O,O)  is  by  definition  equal  to  0.  As  a result,  all  the 
phases  can  be  calculated  from  the  phase  difference  estimates  in 
equation  (1-10). 

Knox  demonstrated  this  technique  with  computer  simulations.  The 
mathematics  of  the  Knox  approach  represent  a dual  analogy  to  that 
underlying  the  motivation  for  the  research  reported  herein. 

Other  approaches  toward  estimating  e^(u,v)  involve  theoretical 


and  in  some  cases  numerical  results  concerned  with  relationships 


1 


between  the  magnitude  and  phase  of  the  OTF  [33-42].  Although  these 
theoretical  results  are  mainly  one  dimensional  and  attempts  to  utilize 
the  theoretical  results  in  practical  algorithms  have  concentrated  only 
on  the  most  rudimentary  functions,  several  examples  will  be  mentioned 
to  indicate  the  present  state  of  these  approaches  which  attempt  to 
extract  phase  information  from  magnitude  information. 

For  example,  Walther  [34]  proves  the  following  theorem. 

Let  h(^}  be  a square  integrable  function  in  the  interval 
and  let  it  be  zero  outside  of  this  interval.  Let  z = x + jy  be  a 
complex  variable,  and  let 

H(z)  » f**  h(c)e'^^’'^^dc 

-»s 

Let  the  zeros  of  H(z)  be  denoted  by  z^. . Then  any  solution  of  the 
equation 

|H(z)|  = |H(z)| 


for  which  the  unknown  function  H(z)  is  band-limited,  must  be  of  the  form 

jc,+jc,z  1-z/zJ 

H(z)  = e ' ^ H(z)  n 


in  which  c^  and  C2  are  real  constants  and  in  which  the  product  is 

extended  over  arbitrarily  many  zeros  of  H(z). 

Note  that  if  all  of  the  zeros  of  H(z)  are  real,  then  H(z)  is 

jci+jc?z 

unique  up  to  the  phase  factor  e . Examples  of  transfer  functions 

with  real  zeros  include 


Hi(z) 

H2(z) 


= linJii 

z 


« 


z 
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H3  (z)  ■ cos  (z)  . 

Preliminary  results  concerning  attempts  equivalent  to  calculating 
transfer  function  phase  from  transfer  function  magnitude  for  transfer 
functions  with  real  zeros  were  presented  In  reference  [40].  These 
results  are  necessarily  one  dimensional  since  the  theory  has  not  been 
extended  to  two  dimensions. 

Another  approach  Is  presented  In  [39].  Here  Gonsalves  assumes 
knowledge  of  both  |h(x)|  and  |H(u)|  and  demonstrates  two  techniques  on 
one  dimensional  examples.  The  techniques  show  promise;  however, 
theoretical  questions  regarding  uniqueness  and  convergence  are  yet  to 
be  worked  out  nor  have  the  algorithms  been  demonstrated  on  practical 
examples. 

1 .2  Research  Goals 

The  ultimate  goal  of  course  Is  to  take  a blurred  Image  and  with 
only  the  knowledge  of  the  blurred  Image  Itself  remove  the  blur.  In  the 
spatially  Invariant  case  the  realization  of  this  ultimate  goal  has  been 
thwarted  by  lack  of  an  algorithm  which  calculates  the  phase  of  the  OTF, 
0^(u,v).  Presented  In  Chapter  2 Is  a method,  mathematically  dual  to  the 
Ideas  of  Knox  [30-32],  which  shows  promise  of  reaching  forward  closer 
to  this  ultimate  goal.  The  research  herein  Investigates  the  feasibility 
and  extent  to  which  this  method  will  estimate  both  the  phase  of  the 


OTF  and  the  magnitude  of  the  OTF. 


, i 

Chapter  2 

THE  ALGORITHMS  TO  BE  INVESTIGATEO 

The  technique  to  be  studied  attempts  to  remove  degradations  from 
an  Image  using  a minimum  of  knowledge.  Assumed  as  given  will  be  the 
following: 

1)  a blurred  Image, 

2)  the  PSF  Is  spatially  Invariant, 

3)  the  extent  of  the  PSF  Is  small  compared  to  the  extent  of 
the  Image, 

4)  the  Image  Is  not  so  severely  blurred  such  that  one  cannot 
tell  the  general  class,  for  example,  building,  outdoor 
scenes,  etc.,  to  which  the  blurred  Image  belongs. 

5)  and  the  blurred  Image  Is  relatively  noise  free. 

The  emphasis  will  be  on  estimating  the  complex  OTF.  That  Is,  both 

magnitude  and  phase  of  the  OTF.  Once  the  OTF  has  been  estimated,  tech- 
niques known  to  be  successful  given  knowledge  of  the  OTF  will  be  used 
to  estimate  the  undegraded  Image. 

The^general  philosophy  will  be  to  assume  that  all  quantities  are 
continuous,  and  any  discretizations  are  a corruption  of  the  continuous 
process  and  Introduce  errors  Into  the  system.  For  example,  the  Image, 
f(x,y).  Is  assumed  to  represent  a continuous  function  of  Intensity 
while,  likewise,  the  PSF,  h(x,y).  Is  assumed  to  represent  a continuous 
function.  Since  convolutions  of  continuous  functions  are  continuous 
functions,  the  blurred  Image,  g(x,y),  will  also  be  assumed  to  be  a 
continuous  function. 
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Dividing  the  degraded  Image  Into  subimages  which  may  overlap  and 


Indexing  the  subimages  by  1, 


Or  equivalently. 


G^(u,v)  H(u,v)F^(u,v) 


G^(u,v)  ■ H(u,v)F^(u,v)  + E^(u,v) 


(2-1) 


where  E^(u,v)  Is  the  error  Inherent  In  approximation  (2-1). 

Forming  the  product 

G^(u,v)G*(u>au,v>&v)  % H(u,v)H*(u+Au,v+Av)F^(u,v)F*(u+Au,v+Av). (2-2) 
Or  equivalently, 

6^(u,v)6*(u+Au,v+Av)  • H(u,v)H*(u+Au,v+Av)F^ (u,v)F*(u+Au,v+Av) 

+E^(u,V,AU,Av) 

where  E^ (u,v,au,av)  Is  now  used  to  denote  the  error  Inherent  in 
approximation  (2-2). 

One  approach  toward  estimating  H(u,v)  would  be  to  average  over  the 
subimages.  That  1s, 

1 N ^ 

G^(u,v)G^(u+au,v+av)  » H(u,v)H*(u+Au,v+Av) 

1«1 

, N 

• ^ ^ F^(u,v)F^(u+Au,v+av) 

1«1 


* i S E^(u,v,au,av) 
1»1 


(2-3) 


1 ^ * 

Note  that  ^ G^(u,v)G^(u+Au,v+Av)  can  be  calculated  from  the 


I 


degraded  Image.  In  addition  note  that  If  i F4(u,v)F4(u+au,v+av) 
, N N 1 1 

N can  be  estimated,  then  the  products 


H(u,v)H*(u-»-Au,v->-Av)  can  be  estimated.  For  a lossless  degradation, 
which  simply  means  the  volume  under  f(x,y)  Is  equal  to  the  volume 
under  g(x,y), 

H(0,0)  » 1 . 


If  the  products  H(u,v)H*(u+au,v+av)  can  be  estimated  and  given 
H(0,0)  * 1,  rearranging  equation  (2-3)  we  have  a recursive  relation- 
ship where  H*(u+au,v+av)  Is  expressed  In  terms  of  H(u,v): 


H*(u+Au,v+Av) 


1 N , N 

N ]^*3^{u,v)G*(u+Au,v+av)  - ff  E^(u,v,au,Av) 

1=1 


1*1 


1 ^ 

H(u,v)  ^ ^ F^(u,v)F*(u+au,v+av) 
1=1 


(2-4) 


Equation  (2-4)  defines  a recursion  which  could  be  utilized  to  simul- 
taneously estimate  both  magnitude  and  phase  of  H(u,v). 

Another  approach  would  be  to  consider  the  magnitude  and  phase  of 
approximation  (2-2)  separately.  First  considering  the  magnitude 


|G^(u,v)G*(u+Au,v+Av) ( * (H(u,v)H*(u+Au,v+Av) 1 1 F^ (u,v)F*(u+au,v+av) I 


rM 


+ E^(U,V,AU,AV) 


where  E^(u,v,au,av)  denotes  the  error  Inherent  In  taking  the  magnitude 
of  approximation  (2-2).  Analogous  to  equation  (2-4),  we  have 

N , N 


|H*(u+au,v+av) I 


1 _ ^ 1 n 

N S I"  If  E"(u,v,au,av) 

. 1=1  M 


1 N ^ 

lH(u,v)l^^  |F^(u,v)F^(u+Au,v+Av)1  ^2-5) 

1-1 
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1 w ^ 1 M 

Similarly,  assuming  ^ ^ |F^(u,v)F^(u+au,v+av)|  and  ^ e"(u,v,au,Av) 

can  be  estimated,  equation  (2-5)  defines  a recursive  relationship 

whereby  the  magnitude  of  H(u,v)  can  be  calculated. 

Now  considering  the  phases  of  approximation  (2-2),  we  have 

0Q  (u,v)  -0g  (u+AU,V+Av)  = 0^(u,v)  - 0^(u+AU,V+Av) 


+ 0p  (u,v)  - 0p  (u+AU,V+Av) 
''l  '‘i 

-t-  0p  (u,V,AU,Av) 


(2-6) 


where  the  subscripts  6^,  H,  and  F^  correspond  to  their  respective 

functions  and  0p  (u,v,au,av)  denotes  the  error  inherent  in  taking 
“"I 

the  phase  of  approximation  (2-2). 

If  equation  (2-6)  is  averaged  in  some  sense  over  1,  and  if  the 

average  value  of  0p  (u,v)  - 0p  (u+au,v+av)  + 0p  (u,v,au,av)  can  be 
'‘i  M •'i 

estimated,  then  rearranging  equation  (2-6)  and  noting  that  0^(0, 0)  = 0, 
one  may  utilize  equation  (2-7)  to  recursively  estimate  0^(u,v). 

0^(U+AU,V+AV)  = 0^(u,v)  - [0g  (u,v)  - 0g  {U+AU,V+AV^J 


+ [0p  (u,v)  - 0p  (u+AU,V+A^] 
•■i  '‘i 


+ 0p  ru,V,AU,Av; 
*'1 


(2-7) 


where  the  bar  denotes  averaging  in  some  sense  yet  to  be  defined. 

l'  ^ 

In  summary,  1)  assuming  T]  F^ (u,v)F*(u+au,v+av)  and 

1 

u Ej(u,v,au,Av)  can  be  estimated,  the  recursion  defined  by 
^ i-1  ' 

equation  (2-4)  can  be  utilized  to  estimate  the  OTF;  or  2)  assuming 
i ^ 1 F<(u,v)F^(u+au,v+av)  1 and  i ^ E*!'(u,v,au,av)  and 

N i»i  1 1 N 1 


[6c  (u.v)-6t-  (u+Au,v+Av)]  + Op  (u.v.Au.Av)  Can  be  estimated,  equations 
•■i  M ‘'i 

(2-5)  and  (2-7)  can  be  used  to  estimate  the  OTF. 

Note  that  the  first  choice  involves  averaging  complex  numbers 
whereas  the  second  choice  involves  averaging  real  numbers.  On  this 
basis  alone  one  could  conjecture  that  the  second  choice  is  perhaps 
the  preferable  one.  A theoretical  development  by  McGlamery  [43] 
confirms  this  intuitive  notion  that  the  second  choice  is  preferable. 
Further  confirmation  has  been  obtained  experimentally  and  is 
consistent  with  McGlamery' s theoretical  results.  Thus,  hereinafter 
it  will  be  assumed  that  equations  (2-5)  and  (2-7)  define  the  basic 
recursions  to  be  studied. 


i i 
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I Chapter  3 

: NATURE  OF  THE  PHASE 

3.1  Discussion 

[ 

Because  the  emphasis  of  this  research  Is  on  the  study  of  a method 

:!  for  estimating  the  phase  of  the  OTF,  It  Is  perhaps  Instructive  to 

‘/I  ■ > 

I discuss  and  In  some  cases  demonstrate  the  nature  and  properties  of  the 

-jl  phase  In  the  Fourier  transform  domain. 

it  Illustrated  In  Figure  3.1  Is  an  attempt  to  present  a perspective 

U : 

plot  of  the  principal  value  of  the  phase  of  the  Fourier  transform  of  a 
64  X 64  subimage  arbitrarily  chosen  from  a digital  image.  From 

: ( 

j Figure  3.1  It  Is  evident  that  the  phase  corresponding  to  Imagery  Is  a 

I rapidly  changing  and  an  extremely  nonlinear  function. 

! One  property  of  the  phase  of  the  Fourier  transform  which  merits 

i attention  Is  the  Ill-conditioned  nature  of  the  phase.  That  Is,  small 

changes  In  the  spatial  domain  Induce  large  changes  In  the  phase 
domain.  Illustrated  in  Each  of  Figures  3.2a,  3.3a,  3.4a,  and  3.5a 
are  two  functions.  The  first  function  corresponds  to  a 64  pixel  cross- 
section  from  an  arbitrarily  chosen  Image;  the  second  function  is 

^ ; 

Identical  to  the  first  function  except  at  a single  point.  Illustrated 
In  Figures  3.2b,  3.3b,  3.4b,  and  3.5b  are  the  phases  associated  with  the 
Fourier  transforms  of  the  respective  two  functions.  Although  some  of  the 
differences  In  the  two  phase  functions  are  because  of  a possible 
Incorrect  unwrapping,  the  point  to  be  made  Is  fairly  obvious,  especially 

) ‘ 

for  the  higher  frequencies;  the  point  being  that  small  changes  In  a 
function  generally  will  Induce  large  changes  In  the  phase  of  the 
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NORMALIZED  DISTANCE 

(a)  two  functions 


NORMAL IZED  FREQUENCY 


(b)  corresponding  phases 
Figure  3.4  Phases  corresponding  to  two  functions 
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Fourier  transform  of  the  function. 

This  ill-conditioned  behavior  dictates  that  any  technique  of  phase 
estimation  will  be  very  sensitive  to  noise.  In  addition,  possible 
problems  with  the  error  term  in  equation  (2-6),  which  forms  the  basis 
of  the  phase  estimation  technique  investigated  herein,  are  foreshadowed. 

For  the  functions  of  Figures  3.2  through  3.5  the  converse  would  be 
true  also.  Since 

3‘^[3(f)]  = f. 

where  script  3 denotes  the  Fourier  transform  operator,  for  these 
functions  large  differences  in  the  phase  of  the  Fourier  transform 
correspond  to  small  differences  in  the  inverse  Fourier  transform.  With 
this  in  mind  let  us  consider  the  converse;  that  is,  let  us  induce 
errors  in  the  phase  and  plot  the  corresponding  functions.  In 
Figures  3.6  through  3.9  there  are  functions  of  Figures  3.2  through  3.5 
respectively,  together  with  five  other  functions.  The  five  other 
functions  correspond  to  errors  in  the  unwrapped  phase  of  10%,  20%,  30%, 
40%,  and  50%.  Note  the  differences  in  the  plot  scaling  between  Figures 
3.2  - 3.5  and  Figures  3.6  - 3.9.  The  functions  without  phase  errors  are 
plotted  heavier  than  the  other  functions.  For  these  functions  increas- 
ing the  phase  error  increases  the  distor^tion  in  the  functions.  Note, 
however,  the  retention  of  most  of  the  features  of  the  functions  even 
under  severe  phase  distortion. 

It  is  interesting  to  note  that  some  phase  distortions  do  not 
produce  a blurring  effect  in  the  sense  that  one  would  notice  a smearing 
or  fuzzy ing  of  an  image.  An  obvious  example  of  phase  error  producing 
no  smearing  is  where  the  phase  errors  are  the  results  of  shifting  the  ^2 
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Image.  From  elementary  Fourier  analysis  we  have 
» (f(x-a.y-b))  - 

where  f(x-a,y-b)  defines  the  image  f(x,y)  shifted  by  an  amount  (a,b). 

Here  there  is  a phase  distortion  but  no  blurring  in  the  conventional 
sense  of  the  word. 

Let  us  now  consider  an  example  of  a phase  distorted  image. 

Fourier  transforming  the  image  in  Figure  3.10,  multiplying  the  trans- 
form by  the  phase  only  filter  illustrated  in  Figure  3.12,  and  then 
inverse  transforming  the  result.  Figure  3.11  was  obtained.  Although 
the  phase  errors  introduced  were  substantial,  as  evidenced  by 
Figure  3.11,  very  little  smearing  or  blurring  was  obtained.  This  is 
consistent  with  the  previous  results  whereby  a substantial  phase 

i 

distortion  in  the  Fourier  domain  induces  minimal  distortion  in  the  i 

picture  domain. 

The  importance  of  the  phase  information  should  not  be  minimized,  * 

however.  For  example,  if  one  Fourier  transforms  an  image,  retains  the 
magnitude  of  the  transform  but  replaces  the  phases  with  random  phases, 
upon  taking  the  inverse  transform  one  obtains  garbage.  An  example  of 
this  process  is  illustrated  in  Figure  3.13  whereupon  a sharp,  detailed 
image  has  been  reduced  to  nonsense  via  this  process. 

3.2  Phase  Bound 

In  this  section  is  derived  for  the  one  dimensional  case  and  also  ' 

for  the  two  dimensional  case  on  axis  a bound  for  the  phase  of  the 
OTF.  First  let  us  consider  the  one  dimensional  case.  Assume  that  the 
PSF  is  non-negative  and  Is  identically  zero  outside  of  an  Interval  27 


-v  , fm.-p  -.*«  a*vr-.-f» 


[~x^.x^].  Assume  additionally  that  the  PSF  Is  sampled  at  the  points 

N N N ' 

- j j +1 . ,-l ,0,1 -1  and  that  the  samples  -L,-L+l ,. ..,-1,0, 

1,...,L-1,  L correspond  to  the  samples  for  x In  the  Interval  [-x^,x^] 

as  In  Figure  3.14.  Now  consider  the  Fourier  series 


2]  h(k)« 


•(j2Ttku/N) 


h(k)i 

k— L 


-(j2iTku/N) 


^ h(k)cos  -j  ^ h(k)s1n 


= R(u)  - jl(u) 

Note  that  0(u)  “ tan'^  l'®|  • 

One  can  represent  H(u),  for  example,  as  In  Figure  3.15  where  the 
abscissa  corresponds  to  the  axis  associated  with  the  cosine  portion 
of  the  transform,  R(u),  the  ordinate  corresponds  to  the  axis  associated 
with  the  sine  portion  of  the  transform,  I(u),  a unit  of  distance 
along  the  curve  Is  du,  and  to  each  point  on  the  curve  corresponds  a 
value  u together  with  R(u)  and  I(u).  e(u}  Is  the  angle  formed 
between  the  abscissa  and  the  point  on  the  curve  corresponding  to  u. 

Because  the  PSF,  h(x).  Is  nonnegative,  and  ruling  out  the  case 
where  h(x)«0  everywhere,  the  point  on  the  curve  corresponding  to  u>0 
will  be  on  the  positive  R(u)  axis.  Notice  that  for  e(u)  to  be 
greater  than  j or  less  than  - J the  cosine  portion  of  the  transform, 
R(u),  must  go  through  zero.  This  simple  Idea  forms  the  basis  of  the 


following  proof.  For  one  period  of  H(u)  the  maximum  number  of  zeros  of 
R(u)  will  be  derived.  To  each  zero  of  R(u)  a change  of  phase 
of  will  be  assumed.  Thus,  the  phase  cannot  exceed  the  number  of 
zeros  of  R(u)  times  +ir. 

One  anomalous  case  should  be  ruled  out  first.  If  the  transform 
is  as  in  Figure  3.16,  the  above  reasoning  would  not  apply.  That  is, 
if  R(u)  is  zero  over  an  interval,  one  would  obtain  an  infinite  number 
of  zeros.  Let  us  consider  the  transform 


»(“)  ’ i: 


where  u is  a real  variable.  Now  consider  the  transform 


R(2)  . 2 


where  z is  a complex  variable,  i.e.,  z = u+jv.  Note  that  R(z)  is  an 
analytic  function  everywhere.  A well  known  theorem  of  complex  variable 
theory  states  that  assuming  f(z)  is  analytic  in  a region,  unless  f(z) 
is  everywhere  zero  in  the  region,  then  the  zeros  of  f(z)  in  the  region 
are  isolated  [44].  Thus,  the  anomalous  case  represented  in  Figure 
3.16  is  ruled  out. 

Now  let  us  determine  the  maximum  number  of  zeros  for  one  period 
of  R(u).  The  result  makes  use  of  the  result  proven  by  Natanson  [45] 
that  the  trigonometric  polynomial  (3-1) 

T(u)  = ^ a^  cos  (^)  + bk  sin  (^)  (3-1) 


I 


can  have  at  most  2L  real  zeros  in  the  Interval  ^ » ®ven  If  each 

multiple  root  is  counted  the  number  of  times  it  occurs.  For  complete- 
ness the  proof  of  this  result  is  given  in  Appendix  A. 

Now  consider 

R(u)  ■ £ I'k  cos  (&!»)  . 

k’-L 

Since  the  cosine  function  is  an  even  function,  we  have 

R(u)  - ^ a,cos(?S^) 

k=0 

where  ag  * hg,  a^  = h .j+h^,  = h_2+^'2»*  "*®l  ~ ’ 

R(u)  is  now  in  the  form  (3-1)  where  bj^  = 0,  and  thus,  R(u)  can  have  at 
most  2L  real  zeros  in  the  interval 

Again  the  cosine  function  is  an  even  function,  and  since  the  linear 
combination  of  even  functions  is  even,  R(u)  will  be  even.  As  a result, 
if  R{u^)  is  a zero  of  R,  so  is  R(-u^).  Therefore,  R(u)  can  have  at 
roost  L real  zeros  in  the  interval  ^0,  and  at  most  L real  zeros  in 
the  Interval  ^ ,oj. 

The  same  reasoning  may  be  applied  to  I(u),  i.e., 

1(0)  . ^ 
k»0 


where  bg  ■ hg,  b^  * h^-h_^,  b^  = h2-h_2»* • • \ 


I(u)  is  now  in  the  form  (3-1)  where  a|^*0,  and  accordingly,  unless  I(u) 
is  identically  zero,  I(u)  can  have  at  most  2L  real  zeros  in  the 
interval  ^ , ^j.  Since  the  sine  function  is  an  odd  function  and  a 
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linear  combination  of  odd  functions  is  odd. 


I(u)  = -I(-u)  . 


Thus,  if  I{u^)  is  a zero  of  I(u),  so  is  I(-u^)  a zero.  Again  I{u)  can 
have  at  most  L real  zeros  in  the  interval  j^O,  and  at  most  L real 
zeros  in  the  interval  ^ ,oj. 


and  at  most  L real 


zeros  in  the  interval 


For  example,  consider  the  case  where  N=64  and  L=3.  As  u varies 
from  0 to  32,  R(u)  has  at  most  3 real  zeros,  and  unless  I(u)  is 
identically  zero,  I(u)  has  at  most  3 real  zeros.  From  Figure  3.17  note 
that  0(32)  cannot  exceed  +3ff. 

Given  N and  L,  and  assuming  I(u)  is  not  identically  zero. 


-Ln  < 0 


(?)  ^ • 


and  since  0(u)  = -0(-u) 


-Lit  < 9^-  < Lit 


Fbr  the  case  in  which  I(u)  is  identically  zero  we  have 

-Lit  _<  0 < Ln  and 

-Ln  < 0 j)  4 Ln  . 

Note  that  the  equality  can  be  achieved.  For  example,  if 


m . £ h,  cos 


where  h j^  » hj^  * .5  and  hj^  * 0 for  Ikj^L,  then 


R(u)  * cos  and 


/ \ 1 

1 / / • z^ro  of  I(u) 

^ \ 3 zeros 

of  l[\x)  / 

zeros  of  I(u) 

\ 

^2  zeros  of  R(u) 

Figure  3.17  Example  of  phase  bound  for  L * 3 

Thus,  this  Is  the  best  bound  possible 


Now  consider  the  two  dimensional  case,  on  axis.  We  have 


H(u,v) 


H(u,v) 


But  this  we  recognize  as  the  one  dimensional  case.  Thus 


The  result  along  the  v-axis  Is  analogous.  Thus 


T 


r 
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Chapter  4 
RESTORATION 

Once  an  estimate  of  both  magnitude  and  phase  of  the  OTF  has  been 
calculated,  an  estimate  of  the  undegraded  image  will  be  made.  In  the 
Fourier  domain 

G(u,v)  « H(u,v)F(u,v)  . 

One  possible  estimate  of  F(u,v)  would  be 

F(u,v)  = where  H(u,v) 

H(u.v) 

denotes  the  estimate  of  the  OTF.  This  approach  gives  acceptable 
results  in  some  cases,  but  in  other  cases  zeros  or  near  zeros  in 
H(u,v)  may  provide  unacceptable  results. 

As  a general  purpose  filter  which  will  give  generally  acceptable 
results,  the  Wiener  filter  defined  in  equation  (4-1)  is  advantageous. 

H*(u,v)'P-(u,v) 

R(u,v)  = 2 (^-1) 

|H{u,v)l  cp^(u,v)tpp(u,v) 

where  cp^(u,v)  andcp^(u,v)  denote  power  spectral  densities  of  the 
undegraded  image  and  noise  respectively. 

The  Wiener  filter  is  optimal  in  a minimum  mean  squared  error 
sense.  Even  though  the  correlation  between  a minimum  mean  squared 
error  criteria  and  the  subjective  quality  of  imagery  is  not  high, 
the  general  purpose  nature  of  the  Wiener  filter  together  with  the 
acceptability  of  results  obtained  via  Wiener  filtering  resulted  in 
its  exclusive  use  in  this  research. 


38 


IT 


Assuming  the  degraded  Image  g Is  stationary  In  the  wide  sense  [47], 
the  power  spectral  density  of  the  degraded  Image, 05^(0, v).  Is  the 
denominator  In  equation  (4-1),  I.e., 

cpg(u,v)  = |H(u,v)|^3,^(u,v)  +(p^(u,v)  . 

Accordingly,  the  denominator  In  the  Wiener  filter  was  estimated 
from  the  degraded  Image  using  a two  dimensional  extension  of  the 

! 

method  documented  by  Welch  [46].  A prototype  can  be  used  to 

j 

estimate  i’x(u,v).  s 

In  general  a subimage  size  of  64  x 64  pixels  was  used  In  the  cal- 
culations. This  subimage  size  results  In  a Wiener  filter  of  approxi- 
mately 33  X 64  nonredundant  complex  numbers.  Because  the  restorations 
were  to  be  on  Imagery  of  256  x 256  pixels,  and  the  filter  corresponded 
to  Imagery  of  64  x 64  pixels,  several  different  Implementations  of  the 
restoration  filter  were  available. 

One  approach  would  be  to  filter  each  nonredundent  64  x 64  pixel 
Subimage  of  the  degraded  Image. 

Another  approach  would  be  to  obtain  the  convolutional  inverse  of 
the  restoration  filter  and  convolve  the  degraded  Image  with  the  convolu- 
tional Inverse  of  the  restoration  filter.  For  example 

R(u,v)G(u,v)  = R(u,v)H(u,v)F(u,v) 
f(x,y)  = 3‘Mr(u,v)G(u,v)1 

where  f(x,y)  denotes  an  estimate  of  the  undegraded  image  f(x,y).  Thus, 

f(x,y)  = r(x,y)*g(x,y)  . 


This  Implementation  requires  convolving  the  64  x 64  kernal , r(x,y). 
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with  the  degraded  Image  g(x.y). 

Still  another  implementation  would  be  to  interpolate  the  Wiener 
filter  of  approximately  33  x 64  nonredundant  complex  numbers  to  a I 

Wiener  filter  of* approximately  129  x 256  nonredundant  complex  numbers. 

The  interpolated  Wiener  filter  would  then  correspond  to  an  image  of  j 

256  X 256  pixels.  This  implementation  is  diagramned  in  Figure  4.1.  | 

A comparison  of  the  three  implementations  has  been  made.  In  | 

Figure  4.2b  is  the  resulting  image  after  convolving  the  image  in  j 

Figure  4.2a  with  a 5 x 5 matrix  of  elements  of  value  1/25,  i.e.,  a 

lossless  PSF.  Assuming  knowledge  of  the  PSF  and  filtering  each  64  x 64 
subimage  of  the  degraded  image,  the  image  of  Figure  4.2c  was  obtained. 

The  severe  artifacts  at  the  subimage  boundaries  are  apparent. 

The  reason  for  the  artifacts  merits  attention.  When  a discrete 
Fourier  transform  (OFT)  is  calculated,  the  function  being  transformed 
is  assumed  to  be  periodic;  as  a result,  in  all  likelihood  the  periodic 
version  of  the  functionwill  contain  jump  discontinuities  at  the  ends 
of  the  period.  These  jump  discontinuities  induce  artificial  high 
frequencies  in  the  Fourier  domain  and  upon  multiplication  by  the  filter 
and  subsequent  inverse  transforming  cause  severe  artifacts  in  the 
restored  imagery.  These  artifacts  are  usually  located  at  sharp  edges  I 

in  the  imagery  and  appear  as  waves  propagating  away  from  the  edge  ' 

boundary.  Sometimes  this  phenomena  is  referred  to  as  ringing. 

Shown  in  Figure  4.3a  is  the  result  of  using  the  same  filter  used 
to  obtain  the  restoration  in  Figure  4.2  except  the  filter  has  been 
interpolated  to  correspond  to  a 256  x 256  image.  The  ringing  around 
the  image  boundary  is  still  bothersome.  Multiplying  the  degraded 
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Best 

Available 

Copy 


(a)  with  no  rounding  of  the  boundary  edge 


(b)  with  rounding  of  the  boundary  edge 


Figure  4.3  Restoration  using  approach  diagrammed  in  Figure  4.1 
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Image  by  a window  which  rounds  the  outer  boundary  of  the  Image,  results 
in  eliminating  the  ringing  at  the  boundary  of  the  restored  image  at 
the  expense  of  removing  some  of  the  Information  at  the  boundary.  This 
Is  Illustrated  In  Figure  4.3b. 

Comparison  of  the  convolutional  inverse  approach  and  the  approach 
diagrammed  In  Figure  4.1  can  be  seen  In  Figures  4.5  and  4.6.  Figure 
4.5a  Is  a test  pattern.  Figure  4.5b  Is  the  result  of  convolving  the 
PSF  shown  In  Figure  4.4  with  Figure  4.5a.  Figures  4.5c  and  4.5d  are 
the  result  after  applying  the  known  filter  via  a 64  x 64  element 

i 

convolutional  Inverse  and  also  via  the  approach  diagrammed  In  ' 

Figure  4.1  respectively.  The  results  are  essentially  the  same. 

In  Figure  4.6  Is  the  same  comparison  except  a different  PSF  has  ! 

been  used.  In  this  case  the  discrete  version  of  the  PSF  was  a 9 x 9 I 

matrix  of  elements  of  value  1/81.  Figure  4.6a  Is  the  blurred  Image 
and  Figures  4.6b  and  4.6c  are  the  result  after  applying  the  Wiener 
filter  via  a 64  x 64  element  convolution  and  via  the  approach  which 
Interpolates  the  filter  In  the  Fourier  domain  respectively.  Again 
the  results  are  essentially  the  same. 

The  OTF  corresponding  to  the  PSF  diagrammed  In  Figure  4.4  Is 
essentially  nonzero  near  the  frequency  axes.  Near  the  axes  Is  where 
most  of  the  energy  In  the  frequency  domain  of  an  Image  Is  concentrated. 

I 

1 

On  the  other  hand,  the  PSF  used  to  arrive  at  the  degraded  Image  In  I 

1 i 

Figure  4.6a  has  several  zeros  on  axis.  The  absence  of  ringing  in  the  I 

restoration  In  Figure  4.5  together  with  the  presence  of  ringing  In 
the  restoration  In  Figure  4.6  can  be  explained  by  the  absence  of 

1 

zeros  In  the  OTF  assumed  In  Figure  4.5  and  the  presence  of  zeros  In 
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(a)  test  pattern 


(c)  restoration  via  convoluti 
/ a1  inverse 

t 


(b)  degraded  test  pattern 


- (d)  restoration  via  interpolation 
of  filter 


Figure  4.5  Comparison  of  approaches  to  restoration 
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(a)  degraded  test  pattern 


(b)  restoration  via  convolutional  (c)  restoration  via  interpolation 
inverse  of  filter 


Figure  4.6  Comparison  of  approaches  to  restoration 
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Chapter  5 

IMPLEMENTATION  AND  SYSTEM  PERFORMANCE 


5. 1 Introduction 

Equations  (5-1)  and  (5-2)  reiterate  the  two  algorithms  to  be 


analyzed. 


|H*(u+Au,v+Av) I = 


with  H(0,0)  = 1, 


]C  |G^(u,v)G*(u+Au,v+av)I  - ^ E^(u,v,Au,Av)J 


1 ^ 

(u,v)l  ^ 2 IF^(u,v)F*(u+au,v+Av)( 


(5-1) 


0^(U+AU,V+AV)  = 0^(U,V)  - [Og  (u,v)-0Q  (U+AU,V+AV] 
+ [0n  (u,v)  - 0c  (U+AU,V+AV)] 

•■l  M 


+ 0n  (U.V.AU.A'^ 


(5-2) 


with  0^(0, 0)  = 0. 

The  success  or  failure  of  the  above  algorithms  In  estimating  the 
magnitude  and  phase  of  the  OTF  for  real  world  blurred  Imagery  Is 
contingent  on  at  least  five  conditions:  1)  whether  or  not  the  mathe- 
matical model  Is  adequate,  2)  whether  or  not  the  ^nknown  quantities 

1 * '' 

^ |F^(u,v)F^(u+au,v+av) I In  the  magnitude  estimate  and  the  unknown 

quantities  [0p  (u,v)-0p  (u+au,v+av)]  In  the  phase  estimate  can  be 
'"i  M 

adequately  estimated,  3)  whether  or  not  the  averages  adequately 

converge,  4)  whether  or  not  the  Iterations  defined  by  equations  (5-1) 

and  (5-2)  are  stable,  and  5)  whether  or  not  the  error  terms 
, N „ 

liT  Z)  Ev(u,v,au,av)  and  e^  (u,v,Au,Av)  in  equations  (5-1)  and  (5-2) 

" 1-1  ' “"I 


respectively  are  negligible.  Assuming  one  does  not  incur  the  sub- 


stantial expense  of  unwrapping  phases  and  considers  only  the  principal 


value  of  the  phase,  an  additional  contingency  upon  which  the  success  or 


failure  of  the  phase  estimate  depends  is  whether  or  not  one  can  ade- 


quately overcome  the  problem  of  considering  only  the  principal  value  of 


the  phases  0-  (u,v)  and  0^  (u,v). 
F. 


5.2  On  the  Assumed  Model 


For  the  error  free  case  experience  has  shown  that  the  model  of 


yo 

the  degrading  process,  i.e.,  g(x,y)  = jj  h(x,y,e,n)f(e,n)dGdn,  is 
reasonable,  at  least  reasonable  to  the  extent  of  enabling  one  to 


improve  real  world  blurred  imagery.  The  results  of  Cannon  [26-28]  on 


real  world  blurred  imagery  together  with  the  results  to  be  presented  in 


Chapter  7 support  this  conclusion. 


5.3  On  the  Question  of  Estimating  Statistical  Properties  of  the 
Unknowns/Convergence 


To  what  extent  can 


1 ^ _ 
jj  |F^(u,v)F*(u+au,v+av)  1 and  [0p  (u,v)  - 0p  i 

i«l  ^ ^ 


U+AU,V+AV] 


be  estimated? 


Using  prototype  images  which  were  statistically  similar.  Cole  [25] 
and  Cannon  [26]  obtained  estimates  of  ^ 5^  |F^(u,v)l  and  cpp(u,v) 
respectively.  A similar  approach  has  been  studied  herein.  Figure  5.1 


contains  four  statistically  similar  images,  statistically  similar  in  the 


sense  that 


1 N ♦ 

nE  (u,v)F^  (u+au,v+av) 


I 


j 

* ^ 

. ^B. 

mm 

‘ ilfc-  'imm 

I 

is  similar  for  k=l,2,3  and  4 where  the  index  k denotes  each  of  the  four  s 

1 

images.  Note  that  the  subject  content  of  one  image  is  different  from 
the  other  images.  The  difference  in  subject  matter  was  intentional 
to  emphasize  the  point  that  to  different  image  subject  matter  may 
correspond  similar  statistics.  Figures  5.2  and  5.3  contain  the 
results  of  the  corresponding  calculations  of 

1 N 
i=l 

for  the  four  images.  Figure  5.2  assumes  au=1  in  units  of  the  funda- 
mental frequency  and  Av=0.  Figure  5.3  assumes  Au=0  and  Av=l.  Note  the 
similarity  of  the  results  for  the  four  images.  The  results  were  ob- 
tained with  Parzen  windowing  (to  be  explained  later),  subimages  of 
64  X 64  pixels  and  a 50%  overlap  of  subimages.  Because  the  images  each 
were  represented  with  512  x 512  pixels,  N in  this  case  was  225.  The 
smoothness  of  the  results  can  be  taken  as  evidence  that  the  averages 
have  reasonably  converged. 

For  extimating  [6^  {u,v)  - (u+au,v+av)]  prototypes  are  not  of 
any  value.  This  can  be  demonstrated  by  calculating  the  correlation 
coefficient,  r,  between  [0r  {u,v)  - (u+Au,v+av)]  for  two  different 
images  of  the  same  prototype  class.  In  this  case  the  correlation  coeffi- 
cient measures  the  degree  of  tendency  for  small  (large)  average  phase 
differences  corresponding  to  one  prototype  to  correspond  to  small  (large) 
average  phase  differences  corresponding  to  a second  prototype.  The 
correlation  coefficient  ranges  between  -1  and  +1.  If  r=-l,  small 
values  in  one  set  corresponds  to  large  values  in  the  second  set  and 
vice  versa.  If  r*0,  there  is  no  correspondence  between  the  two  sets, 
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(a)  autocorrelation  correspond-  (b)  autocorrelation  correspond 
ing  to  Figure  5.1a  Ing  to  Figure  5.1b 


Figure  5.2  Log,«  of  magnitude  autocorrelation  for  au*1,Av»0 


(c)  autocorrelation  correspond 
Ing  to  Figure  5.1c 


(d)  autocorrelation  correspond 
Ing  to  Figure  5. Id 


) 


li- 


(a)  autocorrelation  correspond- 
ing to  Figure  5.1a 


(b)  autocorrelation  correspond' 
ing  to  Figure  5.1b 


(c)  autocorrelation  correspond- 
ing to  Figure  5.1c 


(d)  autocorrelation  correspond' 
ing  to  Figure  5. Id 


Figure  5.3  l-og,«  of  magnitude  autocorrelation  for  au=0,av=1 


and  if  r=l,  small  values  in  one  set  correspond  to  small  values  in  the 
second  set  and  large  values  in  one  set  correspond  to  large  values  in 
the  other  set. 

Shown  in  Table  5.1  are  the  correlation  coefficients  of  the  average 
phase  differences  for  the  six  pairs  of  the  four  images  in  Figure  5.1. 
There  are  several  cases  to  be  considered.  Assuming  the  subimages  are 
not  multiplied  by  a window  prior  to  the  Fourier  transform  operation, 
the  results  in  Table  5.1a  are  obtained.  Here  the  correlation 
coefficient  is  calculated  from  the  average  phase  differences  in  the 
entire  frequency  domain. 

Again  assuming  the  subimages  are  not  multiplied  by  a window  prior 
to  the  Fourier  transform  operation,  but  including  only  the  on  axis 
average  phase  differences  in  the  calculation  of  the  correlation 
coefficient,  the  results  in  Table  5.1b  are  obtained. 

When  the  subimages  are  multiplied  by  a lag  window  prior  to  the 
Fourier  transform  operation,  the  results  in  Table  5.1c  are  obtained. 

In  this  case  a Parzen  window  was  used,  and  the  correlation  coefficient 
was  calculated  across  the  entire  frequency  domain. 

Note  in  all  cases  the  complete  lack  of  correlation.  To  obtain  the 
average  phase  differences  the  technique  outlined  in  Section  5.6  was 
used. 

These  results  dictate  that  another  approach  is  necessary.  Shown 
in  Figures  5.4  and  5.5  are  the  histograms  of  the  phase  differences 
obtained  when  the  four  images  of  Figure  5.1  are  divided  into  sub- 
images, not  multiplied  by  a window,  and  the  results  Fourier  transformed. 
The  histograms  are  across  the  subimages  and  also  across  the  frequency 
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Table  5.1  Correlation  coefficients  between  average  phase  differences  of 
pairs  of  Images  of  Figures  5.1 


Unwindowed  Case 
Image  1 Image  2 au  av 


1 0 -.06 


1 0 -.04 


1 -.001 


Unwindowed  Case  (on  axis) 
Image  1 Image  2 Au  av  r 


1 0 -.02 


1 -.03 

1 -.02 
1 -.03 


1 .21 
1 .06 

1 -.16 

1 .10 


Windowed  Case 
Image  1 Image  2 Au  Av 


1 0 -.03 


1 0 -.03 

1 0 -.02 

1 0 .04 

0 1 .02 

0 1 -.03 

0 1 .02 

0 1 -.08 

0 1 -.10 

0 1 .08 
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PHfiSE  DIFFERENCE  IN  DEGREES 


(a)  histogram  corresponding  to  Figure  5.1a 
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PHfiSE  DIFFERENCE  IN  DEGREES 

(b)  histogram  corresponding  to  Figure  5.1b 

Figure  5. A Phase  difference  histograms  for  case  of  no 
windowing,  au*1.av*0 
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(c)  histogram  corresponding  to  Figure  5.1c 
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PHASE  DIFFERENCE  IN  DEGREES 

(d)  histogram  corresponding  to  Figure  5. id 

Figure  5.4  (Cont'd)  Phase  difference  histograms  for  case  of  no 

windowing,  au*1.av-0 
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(a)  histogram  corresponding  to  Figure  5.1a 
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PHRSE  DIFFERENCE  IN  DEGREES 

(b)  histogram  corresponding  to  Figure  5.1b 

Figure  5.5  Phase  difference  histograms  for  case  of  no  Mindowing 
au«0,Av»l 
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domain.  The  histograms  In  Figure  5.4  correspond  to  phase  differences 
for  one  step  in  the  u direction;  I.e.,  4u>‘1,Av«0.  Likewise,  the 
histograms  in  Figure  5.5  correspond  to  phase  differences  for  one 
step  In  the  v direction;  I.e.,  au»0,av»1.  The  subimages  were  of 
64  X 64  pixels  and  a 50%  overlap  was  used.  This  Implies  N >>  225. 

It  Is  evident  from  the  histograms  In  Figures  5.4  and  5.5  that  the 
modes  of  the  phase  differences  obtained  are  roughly  180°.  Assuming 
the  histograms  are  normalized  to  unit  area,  the  normalized  histograms 
can  be  taken  as  approximations  to  probability  density  functions  of 
the  phase  differences. 

If  one  assumes  that  when  forming  averages  that  the  samples  form- 
ing the  average  each  come  from  the  same  probability  density  function 
and.  In  addition,  are  uncorrelated,  then  the  mean  of  the  average  Is 
the  same  as  the  mean  of  the  given  probability  density  function  and 
the  variance  of  the  average  Is  related  to  the  variance  of  the 
probability  function  by  the  following  relationship  [47]: 

2 

2 _ °PDF 

°AVE  ■ ~ir 

2 2 
where  Is  the  variance  of  the  average,  Opjjp  Is  the  variance  of  the 

underlying  probability  density  function,  and  N Is  the  number  of 

samples  forming  the  average. 

The  phase  differences  for  a given  frequency  but  corresponding  to 
the  Fourier  transform  of  different  subimages  can  in  all  likelihood  be 
considered  as  Independent  and  at  least  be  considered  uncorrelated.  i 

Thus,  one  may  use  the  histograms  In  Figure  5.4  and  5.5  to  estimate  j 
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convergence  properties  of  the  averages. 

It  Is  apparent  that  the  variances  corresponding  to  the  histograms 
In  Figures  5.4  and  5.5  are  considerable.  As  a result.  In  general, 
for  the  case  where  the  subimages  are  not  multiplied  by  a lag  window 
prior  to  Fourier  transforming,  convergence  of  the  average  phase 
differences  will  be  rather  slow.  As  a practical  matter  this  case  Is 
somewhat  overstated  because  the  histograms  are  across  the  entire 
frequency  plane.  If  one  considers  only  those  frequencies  along  the 
axes,  one  obtains  dramatically  different  histograms.  Shown  In 
Figures  5.6  and  5.7  are  the  histograms  of  phase  differences  corres- 
ponding to  the  Images  In  Figure  5.1  for  only  the  on  axis  frequencies 
and  frequency  steps  of  au*l,  Av*0  and  Au=0,  Av=l  respectively.  Thus, 
for  the  case  where  the  subimages  are  not  multiplied  by  a window  before 
Fourier  transforming,  the  convergence  of  the  phase  difference 
averages  will  be  more  rapid  on  axis. 

Now  consider  the  case  where  the  subimages  are  multiplied  by  a 
window.  In  particular  a Parzen  window,  prior  to  being  Fourier  trans- 
formed. Shown  In  Figures  5.8  and  5.9  are  the  histograms  of  the  phase 
differences  obtained  when  the  four  Images  In  Figures  5.1  are  divided 
Into  subimages, windowed  with  a Parzen  window,  and  the  results  Fourier 
transformed.  The  histograms  are  across  the  subimages  and  also  across 
the  frequency  domain  as  was  the  case  for  the  histograms  In  Figures  5.4 

I 

and  5.5.  The  histograms  In  Figure  5.8  correspond  to  phase  differences 
for  one  step  In  the  u direction,  and  the  histograms  In  Figure  5.9 
correspond  to  phase  differences  for  one  step  In  the  v direction.  Note 
that  the  mode  for  the  windowed  case  Is  0,  whereas  the  mode  Is  180**  for 
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PHASE  DIFFERENCE  IN  DEGREES 

(c)  histogram  corresponding  to  Figure  5.1c 


(d)  histogram  corresponding  to  Figure  5. Id 

Figure  5.6  (Cont'd)  Phase  difference  histograms  for  case  of  no 

windowing,  on  axis,  au*1,  Av*0 
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(a)  histogram  corresponding  to  Figure  5.ia  i 
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(b)  histogram  corresponding  to  Figure  5.1b 

Figure  5.7  Phase  difference  histograms  for  case  of  no  windowing, 
on  axis,Au*0,  av=1 
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PHRSE  DIFFERENCE  IN  DEGREED 


PHASE  DIFFERENCE  IN  DEGREES 

(d)  histogram  corresponding  to  Figure  S.ld 


Figure  5.7  (Cont'd)  Phase  difference  histograms  for  case  of  no 

windowing,  on  axis,  Au-0,  Av*1 
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PHASE  DIFFERENCE  IN  DEGREES 


(a)  histogram  corresponding  to  Figure  5.1a 


(b)  histogram  corresponding  to  Figure  5.1b 

Figure  5.8  Phase  difference  histograms  for  case  of  Parzen  windowing, 
AU"1 , AV*0 
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(c)  histogram  corresponding  to  Figure  5.1c 
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(d)  histogram  corresponding  to  Figure  5. Id 

Figure  5.8  (Cont'd)  Phase  difference  histograms  for  case  of 

Parzen  Mlndowlng,  Au-l,av*0 
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(a)  histogram  corresponding  to  Figure  5.1a 
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(b)  histogram  corresponding  to  Figure  5.1b 

Figure  5.9  Phase  difference  histograms  for  case  of  Parzen  windowing 
au"0,av«1 
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(c)  histogrsm  corresponding  to  Figure  5.1c 
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(d)  histogram  corresponding  to  Figure  5. Id 

Figure  5.9  (Cont'd)  Phase  difference  histograms  for  case  of  Parzen 

wIndOMing,  au*0,Av*1 


the  case  of  no  windowing. 

Listed  In  Table  5.2  are  the  standard  deviations  of  the  histograms 
In  Figures  5. 4-5. 9. 
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5.4a 
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5.6a 

69 

5.8a 
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5.4b 

96 

5.6b 

71 

5.8b 

54 

5.4c 

97 

5.6c 

81 

5.8c 

54 

5.4d 

97 

5.6d 

74 

5.8d 

55 

5.5a 

98 

5.7a 

84 

5.9a 

55 

5.5b 

96 

5.7b 

69 

5.9b 

54 

5.5c 

100 

5.7c 

84 

5.9c 

53 

5.5d 

99 

5.7d 

80 

5.9d 

54 

Table  5.2.  Standard  deviation  In  degrees  of  histograms  In 
Figures  5.4  - 5.9 

Note  that  for  the  Parzen  windowed  cases.  Figures  5.8  and  5.9,  the 
standard  deviations  are  much  smaller;  that  Is,  there  Is  more  correla- 
tion between  adjacent  phases  when  the  subimages  are  multiplied  by  a 
Parzen  window  prior  to  the  Fourier  transform  operation. 

The  average  standard  deviation  for  the  histograms  In  Figures  5.8 
and  5.9  Is  roughly  54®.  Plotted  In  Figure  5.10  Is  the  standard  deviation 
of  an  average  of  N samples  from  an  assumed  underlying  Identically 
distributed  probability  density  function  with  a standard  deviation  of  54°. 
The  standard  deviation  of  the  average  can  be  thought  of  as  the  expected 
error  In  attempting  to  estimate  the  mean  of  the  underlying  probability 
density  function.  For  example.  In  this  case  If  one  desires  to  average 
the  phase  differences  to  the  mean,  and  one  uses  N = 225  froir  Figure  5.10 
< on  the  average  the  estimate  of  the  mean  would  be  off  by  3.5°  . 
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5.4  On  the  Stability  of  the  Iterations 

Iterative  processes  at  times  display  erratic  behavior  known  as 
instability.  The  term  instability  has  somewhat  of  a catchall  nature  in 
that  it  is  used  for  a variety  of  ills.  For  example,  one  possible 
definition  is  where  an  error  at  some  step  in  the  iteration  grows  with 
successive  iterations  to  the  point  at  which  the  error  dominates  the 
result.  Another  example  is  where  an  error  causes  oscillations  about 
the  true  answer. 

Because  the  algorithms  being  studied  are  of  an  iterative  nature, 
stability  considerations  are  of  the  utmost  importance.  Accordingly,  let 
us  consider  in  turn  the  magnitude  recursion  defined  by  equation  (5-1) 
and  the  phase  recursion  defined  by  equation  (5-2). 

Rewriting  the  recursion  in  a simpler  form,  the  recursion  correspond- 
ing to  equation  (5-1)  is  as  follows: 


(5-3) 


where 


= |H(u,v)|  , 

^n+1  * |H*(u+au,v+Av)| , 

^ N * 

AFn  =•  If  2 |F^(u,v)F^(u+Au,v+&v)l, 

i=l 

aGn  = If  |G^(u,v)gJ(u+au,v+av)1 


1 M 

- If  ^ e"(u,v,au,av) 
i=l 


First  let  us  consider  the  effects  of  roundoff  error.  Thus, 
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Intuitively,  one  may  observe  from  equation  (5-4)  that  the  roundoff 
error  will  have  a tendency  to  cancel  out.  More  formally,  taking 
expectation  over  the  probability  density  function  associated  with  the 
roundoff  error  and  assuming  the  roundoff  error  at  step  n Is  uncorrelated 
with  the  roundoff  error  at  step  m,  we  have 
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since  the  mean  of  the  roundoff  error  will  be  zero.  Accordingly,  at 
least  with  respect  to  roundoff  error  no  instabilities  are  present. 

To  estimate  a posteriori  the  magnitude  of  H one  will  not  know  aF^ 
in  the  above  notation  but  must  estimate  aF^  by  using  a prototype  image. 
Thus,  Instead  of  aF^  in  equation  (5-3)  the  iteration  will  be 


^ ' WEK  • 

Let  us  now  assume  an  error  of  10%  in  approximating  aFq  by  aPq, 
say  +10%;  assume  no  other  errors  and  let  us  consider  the  propagation  of 
this  error  on  successive  iterates. 


AGo  aGq 


AGn 

aT^ 


2''  • ' 


= 1.1  H, 


.91  H, 


Note  the  oscillation  of  the  estimated  magnitude  about  the  true 


magnitude.  This  simple  analysis  suggests  a possible  Instability  of  an 
oscillatory  nature. 

Experimental  evidence  confirms  the  possible  instability.  Shown  in 

Figure  5.11a  Is  a perspective  plot  of  the  log^Q  of  the  estimate  of  the 

magnitude  of  the  OTF  corresponding  to  the  PSF  Illustrated  in  Figure  4.4. 

This  estimate  was  calculated  via  the  relationship  cDg(u,v)  « |H(u,v)r 

•'**^(u,v)  {i.e. , Cannon's  approach  [26-28])  while  assuming  knowledge  of 

'P^(u.v).  In  Figure  5.11b  Is  a perspective  plot  of  the  log^Q  of  an 

estimate  of  the  magnitude  of  the  OTF  using  the  recursion  defined  in 

equation  (5-1).  Here  again  the  known  Image  was  used  In  the  calculations. 

That  Is,  the  quantities  i |Fj(u,v)F?(u+Au,v+4v) 1 were  calculated 

N 1,1  1 1 

from  the  known  Image.  This  estimate  is  not  as  smooth  as  the  estimate 
Illustrated  In  Figure  5.11a  but  is  of  roughly  the  same  accuracy.  As 
far  as  restorations  using  the  estimates  represented  in  Figures  5.11a  and 
b are  concerned,  no  perceivable  differences  are  noticeable.  For  ex- 
ample, In  Figures  5.12b  and  5.12c  are  two  restorations  of  the  degraded 
image  in  Figure  5.12a.  Figure  5.12b  assumes  knowledge  of  both  magnitude 
and  phase  of  the  OTF,  and  Figure  5.12c  assumes  knowledge  of  the  phase 
of  the  OTF  and  uses  the  estimate  of  the  magnitude  of  the  OTF  repre- 
sented in  Figure  5.11b  which  assumed  knowledge  of  the  undegraded  image. 
Thus,  assuming  knowledge  of  the  undegraded  Image  f,  it  is  seen  that  the 
Iteration  (5-1)  Is  relatively  stable. 

On  the  other  hand.  If  a prototype  image  is  substituted  for  the 
unknown  Image,  the  magnitude  Iteration  cannot  cope  with  the  errors 
Inherent  In  the  substitution.  In  Figure  5.13a  Is  an  estimate  of  the 


magnitude  of  the  OTF  corresponding  to  the  PSF  Illustrated  In  Figure  4.4. 
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This  estimate  used  the  method  of  Cannon  and  the  average  power  spectrum 
calculated  from  Figures  5.1b,  c and  d as  an  estimate  of  the  power 
spectrum  of  Figure  5.1a  which  was  the  Image  degraded  by  the  PSF  In 
Figure  4.4.  Figure  5.13b  Is  an  attempt  to  present  the  corresponding 
result  using  the  magnitude  recursion  defined  by  equation  (5-1). 

Now  let  us  consider  the  stability  properties  of  the  phase  Iteration 
defined  by  equation  (5-2).  Simplifying  the  equation  for  Illustrative 
purposes,  we  have 
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n n n 

Looking  at  the  first  few  terms  to  obtain  a flavor  of  the  propaga- 
tion of  the  results,  we  have 
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roundoff  errors. 
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since  the  roundoff  error  will  be  zero  mean.  Thus,  roundoff  errors  will 
not  Induce  any  Instabilities  In  the  phase  Iteration. 

In  addition  from  equation  (5-5)  It  Is  seen  that  errors  will  be 
additive.  As  a result,  the  Iteration  Itself  will  be  stable. 

5.5  On  Errors  of  Approximation 

Let  us  first  consider  the  errors  Involved  In  approximating  a 
continuous  process  by  a discrete  process.  Although  the  functions  In 
the  model  of  the  blurring  process  are  best  considered  as  continuous 
functions,  the  digital  computer  Is  discrete  and  at  some  point  the 
processes  must  be  discretized  to  be  compatible  with  digital  processing. 
Attempting  to  estimate  continuous  functions  from  sampled  versions  Is 
not  something  to  be  taken  for  granted. 

Because  a continuous  function  has  been  sampled  at  discrete  points 
to  form  the  digital  Image,  errors  due  to  the  sampling  may  be  Introduced. 
For  example.  If  the  sampling  rate  Is  too  low,  high  frequency  components 
In  the  Fourier  domain  will  Impersonate  low  frequency  components.  This 
situation  Is  referred  to  as  aliasing.  In  addition  to  aliasing  one  must 
be  award  of  the  problem  referred  to  as  leakage.  Due  to  the  necessity 
of  considering  finite  subimages,  frequency  components  corresponding  to 
these  finite  subimages  will  be  spread  over  or  leaked  to  adjacent 
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frequency  components. 


When  one  extracts  a subimage  from  an  Image  and  desires  the 
continuous  Fourier  transform  (CFT)  as  an  economic  necessity,  what  does 
one  in  general  settle  for?  In  general  one  settles  for  the  discrete 
Fourier  transform  (DFT)  as  an  approximation  to  the  CFT.  The  DFT  is 
defined  in  equation  (5-6)  below. 


F(K.L) 


N-1  N-1 


1 I 

7 E E 


[-j2n{IK+JL)/N] 


(5-6) 


1=0  J=0 


where  I,J,K,L  = 0,1,..., N-1  . 

The  DFT  is  calculated  via  the  fast  Fourier  transform  (FFT)  which  is  an 
efficient  algorithm  for  such  calculations. 

The  DFT  can  be  shown  to  be  an  approximation  to  the  CFT.  If  one 
begins  with  the  definition  of  the  CFT,  truncates  the  limits  of  inte- 
igration  to  (>  j ),  samples  the  integrand  at  equally  spaced  points, 
applies  the  trapezoidal  rule  [48]  to  numerically  integrate,  and  ignores 
the  linear  phase  factor  due  to  shifting  the  lower  limit  of  - ^ on  the 
integral  to  a lower  limit  of  0 on  the  summation,  one  obtains  the  DFT. 

Considering  the  one  dimensional  case  for  simplicity,  at  the 

frequency  K » N-1  it  can  be  shown  that  the  factor  in  the 

Integrand  of  the  definition  of  CFT  is  sampled  at  two  points  per  cycle 

where  cycle  refers  to  the  periodic  nature  of  i.e.,  = 

cos  2irux  - J sin  2wux.  Thus,  at  the  frequency  K ■ N-1  the  function 

to  be  transformed  is  multiplied  by  a function  which  is  changing  rapidly 

relative  to  the  sampling.  Thus,  it  is  not  hard  to  see  that  in  using 

the  FFT  to  approximate  continuous  Fourier  transforms  a degree  of  error 

82 


Is  Introduced. 


Both  f(I>d)  and  F(K»L)  are  now  assumed  to  be  periodic,  and  as  a | 

consequence  of  the  periodicity,  typically  the  function  f(I,J),  now 
considered  as  a periodic  function,  contains  Jump  discontinuities  at  the 
ends  of  the  periods.  The  jump  discontinuities  Induce  spurious  high 

I 

frequency  components  In  the  Fourier  domain.  To  alleviate  this  | 

situation  the  functfon  f(I,J)  Is  multiplied  by  a function  w(I,J),  | 

! 

called  a window,  which  smoothes  to  zero  the  ends  of  the  periodic  | 

function,  thus  eliminating  the  jump  discontinuities  and  the  associated 
sharp  edges.  As  a consequence  of  the  fact  that  multiplication  In  the 
picture  domain  corresponds  to  convolution  In  the  frequency  domain,  we 
would  have  In  the  frequency  domain  upon  multiplying  f(I,J)  by  w(I,J) 

F(K,L)*W(K,L) 

where  * denotes  convolution. 

At  this  point  It  Is  appropriate  to  consider  the  choice  of  the 
window  w(I,J)  for  the  problem  at  hand.  Ultimately,  an  estimate  of 
H(u,v)  Is  desired.  Since  the  blurred  Image  Is  assumed  to  be  a 
convolution  across  the  entire  Image,  we  have  only  the  approximate 
relationship  i 

g^(x,y)  h(x,y)*f^(x,y)  . (5-7)  ; j 

Referring  to  the  width  of  h(x,y)  as  the  portion  of  the  Independent 

j 

variables  corresponding  to  nonzero  h(x,y),  then  for  points  of  g^(x,y)  ' * 

within  a half  width  of  h(x,y)  from  the  boundary  of  g^(x,y),  approxi- 
mation (5-7)  Is  exact.  For  points  of  g^(x,y)  outside  a half  width  of  | 

h(x,y)  from  the  boundary  the  relationship  Is  erroneous.  Note  by  i 
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equation  (1-2)  approximation  (5-7)  assumes  by  definition  an  aperiodic 
convolution.  On  the  other  hand.  In  the  discretized  domain  In  which  by 
necessity  we  must  do  our  calculations,  we  have  the  pair 

g(I,J)  = h(I,J)©f(I,J) 
and 

G(K,L)  = H(K,L)F(K,L) 

where  © denotes  a circular  or  periodic  convolution.  Here  h(I,J)  and 
f(I,J)  are  assumed  to  be  periodic,  and  as  a result,  the  circular  convolu- 
tion can  be  used  as  an  approximation  to  the  aperiodic  convolution  only 
with  caution. 

The  relationship  between  the  aperiodic  and  circular  convolution  Is 
again  for  points  of  g^(x,y)  within  a half  width  of  h(x,y)  from  the 
boundary,  the  two  are  the  same.  For  points  of  g^(x,y)  outside  a half 
width  of  h(x,y)  from  the  boundary,  the  difference  between  the  aperiodic 
and  circular  convolution  can  be  substantial. 

The  differences  between  the  aperiodic  and  circular  convolution  are 
not  only  relevant  as  an  approximation  which  Induces  error  Into  the  cal- 
culations but  additionally  because  the  errors  In  approximation  (5-7) 
occur  at  the  same  points  as  the  errors  Induced  by  replacing  the 
aperiodic  convolution  with  the  circular  ronvolutlon. 

With  this  In  mind  let  us  consider  several  candidate  windows. 
Including  the  case  of  no  window,  and  compare  the  error  In  the  Fourier 
domain  for  functions  convolved  aperlodically  and  windowed  by  a candidate 
window  and  also  convolved  circularly  and  windowed  by  the  same  candidate 


window. 

First  let  us  consider  the  case  of  no  windowing  on  several  one 
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dimensional  cross  sections  picked  arbitrarily  from  an  Image.  Each  of 
these  functions  has  been  convolved  aperlodically  and  circularly  and 
displayed  In  Figures  S.14a»  5. 15a.  5.16a  and  5.17a  respectively.  Note 
from  these  figures  that  the  width  of  the  PSF  Is  7 pixels.  In  Figures 
5.14b,  5.15b.  5.16b  and  5.17b  are  the  phase  functions  In  the  Fourier 
domain  of  the  respective  two  functions  displayed  In  Figures  5.14a,  5.15a, 
5.16a  arid  5.17a.  It  Is  fairly  obvious  that  for  the  case  of  no  windowing 
the  error  term  In  equation  ( 2 >6 ), repeated  here  for  convenience,  can  be 
substantial. 


0-  (u,v)  - 0-  (u+Au,v+Av) 

bi 


= 6j^(u,v)  - 0j^(u+AU,V+Av) 

+ 0r  (U,V)  - 0c  (U+AU,V+AV) 

■f  Oc  (U.V.AU.AV) 

*^1 


This  Is  another  demonstration  of  the  Ill-conditioned  nature  of  the 
phase  of  the  Fourier  transform. 

Next  consider  the  candidate  windows  Illustrated  In  one  dimension 
In  Figure  5.18.  These  windows  are  usually  known  as  Hamming,  Hanning, 
triangle,  and  Parzen  windows.  In  Figure  5.19  Is  a perspective  plot  of 
an  arbitrarily  chosen  64  x 64  pixel  subimage.  Figure  5.20  contains 
ijperspecti ve  plots  of  the  log^Q  of  the  absolute  value  of  the  difference 
In  magnitude  In  the  Fourier  domain  between  the  subimage  illustrated  In 
Figure  5.19  convolved  aperlodically  and  circularly  by  the  PSF  displayed 
In  Figure  4.4  and  windowed  by  the  Indicated  windows.  The  differences 
have  been  normalized  by  the  DC  terms,  (0,0)  frequency  component,  of  the 
aperlodically  convolved  result.  Note  that  for  all  four  windows  the 
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(a)  a function  convolved  aperlodically 
and  circularly 
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(b)  corresponding  phases 

Figure  5.14  Phases  corresponding  to’ a function  convolved  aperlodically 
and  circularly 


(a)  a function  convolved  aperiodically 

and  circularly 
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(b)  corresponding  phases 
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Figure  5.16 


Phases  corresponding  to  a function  convolved  aperiodically 
and  circularly 
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(a)  results  for  triangle 
windOMlng  (highest 
contour  level  = -3) 


results  for  Hanning 
windowing  (highest 
contour  level  *-4) 


(b)  results  for  Hanning 
windowing  (highest 
contour  level  = -3) 


(d)  results  for  Parzen 
windowing  (highest 
contour  level  =-5) 


Figure  5.20  Magnitude  differences  of  Fourier  transforms  of 
subimage  convolved  aperiodically  and  circularly 
for  different  windows  (difference  in  contour 
levels  * 1) 


differences  are  comparable  and  suggest  that  the  approximation  (used  by 
both  Cole  [25]  and  Cannon  [26-28]) 

|G^(u,v)|  % |H(u,v)l|F^(u.v)| 

Is  well  behaved  and  realistic. 

Figure  5.21  contains  the  corresponding  results  for  the  differences 
Mn  phases  In  the  Fourier  domain.  All  four  perspective  plots  have  the 
same  plot  scaling;  thus,  are  directly  comparable.  Ng  attempt  has  been 
made  to  unwrap  the  phases,  and  the  differences  are  actually  the 
absolute  value  of  the  principal  value  of  the  difference  in  principal 
values  of  the  two  sets  of  phases  for  a given  windowing.  Accordingly, 
the  maximum  error  will  be  180®.  Although  the  phases  have  not  been 
unwrapped,  the  point  to  be  made,  however,  is  clear.  Because  the  phases 
in  the  Fourier  domain  are  susceptible  to  large  changes  from  small 
changes  in  the  picture  domain,  I.e.,  ill-conditioning,  the  approximation 

0G^(u,v)  ^ ej^(u,v)  6p^(u,v)  (5-8) 

Is  not  a particularly  good  one  In  general.  Nevertheless,  In  particular 
the  Parzen  windowing  gives  resulting  differences  In  phases  that  are 
relatively  small.  This  Is  due  to  the  fact  that  the  "tails"  of  the 
Parzen  window  are  essentially  zero  over  the  region  In  which  the 
aperlodically  convolved  result  and  circularly  convolved  result  differ. 
Thus,  the  Parzen  windowing  In  effect  throws  away  the  values  In  the 
picture  domain  that  differ  as  a result  of  the  differences  between  the 
aperiodic  and  circular  convolution. 

Since  approximation  (5-8)  Is  a crucial  assumption  leading  to  the 
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(a)  results  for  triangle 
windowing 


(c)  results  for  Hanning 
windowing 


(b)  results  for  Hamming 
windowing 


(d)  results  for  Parzen 
windowing 


Figure  5.21  Phase  differences  of  Fourier  transforms  of  subimage 
convolved  aperiodically  and  circularly  for  different 
windowing. 
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recursive  estimate  of  phase,  that  Is,  equation  (5-2),  and  because  the 
\approx1mat1on  Is  not  a particularly  sharp  approximation,  the  technique 
must  be  reformulated  to  take  advantage  of  the  advantageous  properties 
of  the  Parzen  windowing.  Accordingly,  Including  the  effects  of  the 
window,  we  now  have 

6(u,v)*W(u,v)  = [H(u,v)F(u,v)]*W(u,v) 

where  W refers  to  the  Fourier  transform  of  the  window.  For  a subimage 
we  have  approximately 

G^(u,v)*W(u,v)  % [H(u,v)F^(u,v)]*W(u,v)  (5-9) 

and  the  further  approximation 

G^(u,v)*W(u,v)  ^ H(u,v)[F^(u,v)*W(u,v)]  (5-10) 

can  be  made.  This  approximation  can  be  reasoned  on  the  Idea  that  H(u,v) 
Is  essentially  constant  across  the  extent  of  the  spectral  window  W(u,v) 
and  therefore  can  be  factored  outside  the  Integral  Implicit  In  the 
convolution  operator.  If  one  Is  concerned  with  magnitudes  at  this 
point,  then 

|G^(u,v)*W(u,v) I % |H(u,v) I |F^(u,v)*W(u,v) |.  (5-11) 

Experimental  evidence  confirms  the  appropriateness  of  approximation 
(5-11),  and  as  judged  by  the  results  of  both  Cole  [25]  and  Cannon  [26- 
28]  who  used  this  approximation  In  their  work,  the  approximation  Is 
Indeed  realistic.  As  a result,  the  recursive  estimate  of  OTF  magnitude 
with  the  Inclusion  of  the  window  Is  not  substantially  different  than 
equation  (5-1).  For  completeness  the  analogous  recursion  Is  Included 
as  defined  by  equation  (5-12)  below. 

95 


I 
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|H(u,v)|  TT  £ |[F.(u,v)*W(u,v)][Ft(u+Au,v+Av)*W(u+Au,v+Av)]| 

i=l  1 1 


rM, 


(5-12) 

where  the  error  term  E”(u,v,Au,av)  Is  denoted  by  the  same  notation  as 
in  equation  (5-1)  but  is  in  general  different. 

Let  us  now  include  the  effects  of  the  window  in  the  formulation  of 
the  phase  estimate.  Analogous  to  approximation  (5-11)  we  have 


3g^*y^(u,v)  ^ e^(u,v)  + 0f^*h(u.v) 


where  0g  *y(utv)  denotes  the  phase  of  G^(u,v)*W(u,v)  and  *y(u,v) 
denotes  the  phase  of  F^(u,v)*W(u,v).  This  approximation  leads  to  the 
phase  recursion  (5-13)  . 


0j^(u+Au,v+Av)  = 0j^(u,v)  - 


+ - ®F,  ^y(u+AU,V+Avl]  + E(U,V,AU,AV) 

(5-13) 

where  E denotes  the  error  inherent  in  the  assumed  approximations. 

The  operation  of  convolution  in  the  transfonn  domain  results  in  a 
smearing  or  blurring  of  the  Fourier  coefficients.  One  might  ask  to 
what  extent  can  one  utilize  the  phase  differences  of  the  smeared 
quantities,  i.e.,  0q  *y(u,v)  and  Op  *y(u,v),  in  estimating  the  phase 
differences  of  the  OTF.  To  answer  this  question  the  following 


simulation  was  performed. 
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A 512  X 512  pixel  Image  was  subdivided  Into  225  64  x 64  pixel 
subimages  where  a 50%  overlapping  of  subimages  was  allowed.  Each  of 
the  225  subimages  were  Fourier  transformed  and  the  Fourier  coefficients 
on  or  within  two  fundamental  frequencies  of  the  axis  were  saved  on  tape 
for  later  processing. 

In  the  simulation  the  saved  coefficients  correspond  to  F^(I,J) 
where  I and  J are  associated  with  those  Indices  on  or  within  two  steps 
of  the  frequency  axes.  At  this  point  one  can  assume  a given  OTF, 
multiply  the  given  OTF  times  each  F^.(I,J)  and  obtain 

G^(I.J)  = H(I,J)F^{I.J)  ./  (5-14) 

Note  that  by  definition  relationship  (5-14)  will  be  exact,  and 
computer  Implementation  of  the  multiplication  will  preserve  the  equality 
to  within  the  roundoff  error  qf  the  machine.  Thus,  using  quantities 
calculated  via  equation  (5-14),  one  can  study  the  effects  of  various 
windows  by  the  calculation  of 

G.(I,J)*W(I,J)  . 

The  simulation  assumed  that  the  Fourier  coefficients  along  the 
u-axis  and  the  Fourier  coefficients  along  the  v-axis  had  the  same 
underlying  probability  density  functions,  and  thus  were  lumped  together. 
Accordingly,  the  results  will  be  considered  to  be  along  an  axis.  For 
definiteness  assume  the  I axis. 

Letting  the  phase  of  H be  linear  and  one  dimensional, 

e„(I,J)  = al,  (5-15) 


one  can  calculate 
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In  the  Ideal  case  one  would  desire  the  average  at  each  frequency  I In 
(5-16)  to  be 

al  - a(I-*-l ) = -a. 

Even  if  the  Ideal  was  not  found,  but  If  to  each  "a,"  the 
average  (5-16)  converged  to  a unique  "b,"  then  given  a "b"  the  correct 
phase  difference  of  the  OTF  could  be  found.  That  Is,  If 

b(I)  = f(a(I)) 

where  f denotes  a function  and  assuming  to  each  a(I)  corresponds  a 
unique  b(I),  then 

a(I)  - f^(b(I)) 

where  f”^  denotes  the  Inverse  mapping. 

Shown  In  Figures  5.22  and  5.23  are  the  results  of  the  simulation 

'I  for  the  case  of  no  windowing.  The  subimages  used  In  the  simulation 

f!: 

; ii  corresponding  to  Figure  5.22  were  taken  from  the  Image  shown  In 

1,1  Figure  5.1a;  the  subimages  used  In  the  simulation  corresponding  to 

^ f 

Figure  5.23  were  taken  from  the  Image  shown  In  Figure  5.1b.  In  this 

simulation  e^(I,J)  assumed  nine  different  linear  functions.  The  first 

function  was  (I,J)  - 0.  The  second  function  was  linear  such  that 
”l 

sj  Ou  (I,J)  - 0u  (1+1,0)  * -11.25° 

■ I ^2 

1 

The  third  function  was  linear  such  that 


eu  (1,0)  - Ou  (1+1,0)  * -22.50° 

"3  ”3 

i 

* 

I and  so  on 
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NORMRLIZED  FREQUENCY 

(b)  linear  least  squares  fits  to  functions  in  (a) 

Figure  5.22  Sinxilation  results  corresponding  to  no  windowing  and 
Figure  5.1a 


NORMALIZED  FREQUENCY 


(a)  phase  averages 


NORMAL IZED  FREQUENCY 


(b)  linear  least  squares  fits  to  functions  in  (a) 

Figure  5.23  Simulation  results  corresponding  to  no  windowing  and 
Figure  5.1b 


-90.00° 


r 


Ou  (I.J)  - e.  (I+l.J)  = 

Hg  ng 


Each  curve  in  Figures  5.22a  and  5.23a  corresponds  to  the  average 


e.  (I.O)  - e.  (1+1,0) 

bi  b^ 

corresponding  to  each  of  the  nine  0m  (I)  functions.  Here  i ranged  over 

”k 

225  subimages  and  the  two  axes  which  were  lumped  together.  Each  curve  . 

i 

in  Figures  5.22b  and  5.23b  is  a linear  least  squares  fit  to  the  curves  ! 

'I 

in  Figures  5.22a  and  5.23a  respectively.  Note  from  the  least  squares  fits  i 

that  the  convergence  is  roughly  to  the  correct  phase  differences.  | 

Recall,  however,  that  relationship  (5-14)  has  been  induced  artificially  | 

and  the  problem  inherent  in  approximation  (5-8)  has  not  been  overcome. 

Assuming  the  same  nine  phase  functions  and  assuming  Parzen 
windowing,  the  curves  in  Figures  5.24  and  5.25  were  calculated.  Here 
each  curve  in  Figures  5.24a  and  5.25a  correspond  to  the  average 

corresponding  to  each  of  the  nine  Sm  (I)  functions.  Again  i ranged  over 
the  225  subimages  of  Figures  5.1a  and  5.1b  respectively  and  ranged  over 
the  two  axes  which  were  lumped  together.  Each  curve  in  Figures  5.24b 
and  5.25b  is  a linear  least  squares  fit  to  the  curves  in  Figures  5.24a 
and  5.25a  respectively. 

From  Figures  5.22  and  5.23  and  from  Figures  5.24  and  5.25  it  is 
seen  that  the  results  corresponding  to  the  different  images  of  Figures  ■ 

5.1a  and  5.1b  are  roughly  the  same.  ' 

An  important  observation  drawn  from  Figures  5.24  and  5.25  is  when 
assuming  a Parzen  window,  the  inability  of  the  averages  (5-17)  to 
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(b)  linear  least  squares  fits  to  functions  in  (a) 

Figure  5.24  Simulation  results  corresponding  to  Parzen  windowing 
and  Figure  5.1a 


resolve  phase  differences  of  the  OTF  of  value  less  than  approximately 
40“ . 

In  addition  note  the  phase  differences  in  Figures  5.24  and  5.25 
have  the  incorrect  sign.  With  this  in  mind  let  us  focus  our  attention 
on  the  mechanics  of  forming  a discrete  convolution.  The  definition  of 
discrete  convolution  is  given  in  equation  p-18)  below 

, N-1  N-1 

G{I,J)*W(I,J)  = T S 51  6(I-K,J-L)W(K,L)  (5-18) 

^ K=0  L=0 

For  simplicity  assume  the  width  of  W is  five  units  in  each  direc- 
tion. That  is,  assume  in  a square  of  five  values  by  five  values  the 
values  of  W are  nonzero  and  the  other  values  of  W are  zero.  For 
definiteness  assume  additionally  that  the  nonzero  values  of  W correspond 
to  the  K and  L indices  of  0,1 ,2,N-2,and  N-1.  But  since  periodicity  is 
assumed  this  is  equivalent  to  assuming  the  nonzero  values  of  W corres- 
pond to  theKand  L indices  of  -2,-1 ,0,1,  and  2.  Except  in  unusual 
circumstances  will  W be  anything  but  real  and  symmetric;  accordingly, 
assume  that  W is  real  and  symmetric. 

Under  these  assumptions  let  us  look  at  the  quantity 

[G(I,J)*W(I,J)][G(I+1,J)*W(I+1,J)].*  (5-19) 


Expanding,  we  have 

, 2 2 


(5-20) 


and 


[G(I*1,J)*W(U1,J)]* 


1 ^ 2 

4 V V h*(i*-i-k,j-l)f*(i+i-k.j-l)w(k,l) 

N K-2  L -2 

^ I ^ 104 


T 


Multiplying  the  right  hand  side  of  equation  (5-20)  by  the  right 
hand  side  of  equation  (5-21)  results  in  the  sum  of  625  terms.  Let 
us  group  each  of  the  625  terms  into  two  categories.  The  first  category 
consists  of  terms  which  contain  factors  of  the  form 


H(I.J)H*(I.J)F(IJ)F*(I.J) 

and  the  second  category  consists  of  all  other  terms.  In  this  example 
20  of  the  625  terms  would  be  in  the  first  category. 

The  terms  in  the  first  category  contain  no  phase  information, 
whereas  the  terms  in  the  second  category  contain  the  phase  information 
of  the  product. 

Multiplying  equation  (5-20)  times  equation  (5-21),  we  have 

[G(I,J)*W(I,J)][G(I+1,J)*W(I+1,J)]*  = A(I,J)  + B(I,J) 

where  A is  the  sum  of  those  20  terms  of  the  first  category  and  B is  the 
sum  of  the  remaining  605  terms  of  the  second  category.  Figure  5.26 
illustrates  the  vector  sum  A+B. 

Assuming  the  degraded  image  is  divided  into  subimages,  the  corres- 
ponding results  would  be 

[G^(I.J)*W(I,J)][G^(I+1,J)*W(I+1.J)]*  = A.(I,J)  + B.(I,J)  . 

Recall  from  Figures  5.4-5. 7 that  the  phase  differences  corresponding 
to  an  unwindowed,  undegraded  image  will  on  the  average  be  out  of  phase 
for  one  step  in  the  frequency  domain.  Additionally,  the  terms  implicit 
in  vector  B^  which  predominate  will  be  those  terms  associated  with  one 
step  in  the  frequency  domain.  As  a result,  on  the  average  vector  B. 
will  be  out  of  phase  by  +n. 
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Surprisingly,  on  the  average  the  magnitude  of  A^(1,J)  Is  roughly 
the  same  as  the  magnitude  of  Experimental  evidence  of  the 

calculated  average, 

, 100  |B^(I,0) 

[ 1-1  ’ 

• for  I corresponding  to  the  data  on  axis,  using  Parzen  windowing,  and 

I 100  subimages  of  Figure  5.1a  Is  given  In  Table  5.3. 

k 

I This  explains  the  results  of  Figures  5.4-5. 7 compared  to  the 

[ results  of  Figures  5.8  and  5.9  In  the  sense  that  the  phase  differences 

[ corresponding  to  unwindowed  subimages  tended  to  be  out  of  phase  while 

ji  the  phase  differences  corresponding  to  the  Parzen  windowed  subimages 

i 

|i  tended  to  be  In  phase.  For  example,  assuming  vectors  A and  B are 

roughly  the  same  magnitude,  the  vector  A+B  In  Figure  5.26  will  tend  to 
have  a mean  of  zero  phase  even  though  the  vector  B will  have  a mean  of 
phase. 

In  addition  this  explains  the  sign  discrepancy  In  the  results  of 
Figures  5.24  and  5.25  where  the  Parzen  window  was  used.  Assuming  that 
the  phase  between  Index  I and  Index  I+l  of  the  OTF  Is  Increasing,  the 
phase  difference  o^( I )-o^( 1*1 ) will  be  negative.  The  average  of  the 
vector  B corresponding  to  this  Index  will  be  as  Illustrated  In 
Figure  5.27.  But  the  averages  corresponding  to  Figures  5.24  and  5.25 
are  the  averages 

In  Figure  5.20  Is  an  Illustration  of  the  dynamics  behind  the  results  of 

Figures  5.24  and  5.25.  This  Illustration  assumes  seven  subimages  for 
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• (DC) 

1 

.01 

2 

3 

3.30 

4 

5 

.93 

6 

7 

.99 

8 

9 

.91 

10 

11 

1.01 

12 

13 

1.03 

14 

15 

1.00 

16 

17 

.92 

18 

19 

.98 

20 

21 

.96 

22 

23 

1.02 

24 

25 

.93 

26 

27 

.93 

28 

29 

.88 

30 

31 

.90 

Table  5.3  Average  magnitude  ratio  of  category  two  vectors 
divided  by  category  one  vectors 
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Figure  5.26  Vector  sum  A+B 


Figure  5.27  Average  vector  B corresponding  to  H(I) 


simplicity.  Assume  that  B -is  roughly  the  vector  denoted  by  Bg.  Note 
that  each  of  the  have  zero  phase  and  that  e^(I)-e^(I+l ) is  roughly 
-30°.  The  quantity 


will  be  the  average  phase  of  the  resultant  vectors  A^.+B^. . Note  that 
although  the  vectors  B^ ,62.82, B^.Bg,  and  B^  are  more  or  less  symmetri- 
cally distributed  about  the  vector  Bg,  five  of  the  resultant  vectors 
A^+B^  have  positive  phase  and  two  of  the  resultant  vectors  A^+B^. 
have  negative  phase.  Thus, 


will  be  of  different  sign  than  9^(I)-ej^(I+l ) . Additionally,  from  the 
experimental  results  the  magnitude  of  the  phase  difference  average 
corresponding  to  Parzen  windowed  subimages  will  be  less  than  the 
magnitude  of  the  phase  difference  of  the  underlying  OTF. 

5.6  On  Phase  Averaging 

Forming  a phase  average  presents  difficulties.  Usually  one  only 
knows  the  principal  value  of  the  phase.  Illustrated  in  Figure  5.29  is 
representation  of  the  complex  value  for  an  example  H(u,v),  i.e., 
H(u^,v^).  Here  x and  y refer  to  the  complex  number  z=x+jy.  Normally 
one  does  not  know  whether  the  phase  of  H(upV|)  is  9^(U|,v^)  as 
illustrated  or  9^(u^,Vi)  + 2iri  where  i is  an  integer. 

Using  the  definition  (5-22)  as  the  definition  of  the  average  may 
lead  to  incorrect  results. 


no 


» = R i:»t  (5-22) 

1-1 

For  example,  the  phases  of  the  vectors  in  the  complex  plane  of  Figure 
5.30  obviously  average  to  around  n.  Assuming  a principal  value  of  0 to 
2n  in  definition  (5-22),  one  obtains  the  correct  average  of  n.  On  the 
I other  hand,  if  one  assumes  a principal  value  of  -n  to  n in  definition 

(5-22),  then  the  average  phase  will  be  an  erroneous  0. 

Because  the  Fourier  transform  of  a continuous  function  is  continu- 
ous, the  real  and  imaginary  parts  of  the  Fourier  transform  are  continu- 
ous. As  a result,  for  the  most  part  the  phase  of  the  Fourier  transform 
is  continuous.  The  exception  being  possible  jump  discontinuities  of 
_+7i  when  the  complex  valued  Fourier  transform  goes  through  0. 

"Unwrapping"  the  phase  refers  to  determining  the  continuous  phase  where 
the  possible  jump  discontinuities  of  +it  are  allowed.  If  one  assumes 
unwrapped  phases,  then  the  average  in  definition  (5-22)  is  well  defined. 

Techniques  for  estimating  the  unwrapped  phase  have  been  studied 
by  other  researchers  [29,49].  These  techniques  necessitate  increased 
computation  and  an  alternative  approach  was  taken. 

Referring  to  equation  (2-6)  and  neglecting  the  error  term,  we  have 


(u,v)  - 0- 

Gi 


(u+AU,V+Av)  ^ o^(u,v)  - e^(u+AU,V+Av) 


+ 3r  (u,v)  - 0p  (u+AU,V+Av) 


The  phase  difference  corresponding  to  the  degraded  subimage  is  the  sum 
of  the  phase  difference  of  the  OTF  and  the  phase  difference  correspond- 
ing to  the  undegraded  subimage.  For  a fixed  u,v,Au,  and  Av,  say  u*Ui , 
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v-v^,  AU»AU^,  and  av»av^,  e^(u^,v^)  - e|^(u^+Au^ ,v^+av^ ) will  be  fixed 
for  each  subimage  and 


Op  (u^.v^)  - Op  (u^+Au^,v^+Av^) 


1 ■ ■ '1 

will  vary  from  subimage  to  subimage.  For  example,  In  Figure  5.31a  Is 
represented  a unit  vector  of  phase  equal  to  e^(u^,v^)  - 0^(u^+AUpV^+AV|); 
In  Figure  5.31b  Is  represented  unit  vectors  corresponding  to 
0p  (u^,v^)  - 0p  (u^  +Au^,v^+Av^)  associated  with  several  subimages; 
Figure  5.31c  represents  the  unit  vectors  obtained  by  summing  the  phase 
difference  In  Figure  5.31a  to  each  of  the  phase  differences  in 
Figure  5.31b.  Note  that  the  vectors  In  Figure  5.31c  are  the  vectors 
In  Figure  5.31b  rotated  by  the  phase  of  the  vector  in  Figure  5.31a. 

If  the  phase  differences  In  Figure  5.31b  are  of  zero  mean,  then 
the  mean  of  the  phase  differences  In  Figure  5.31c  should  be 


0^(u.|,v^)  - 0^(U^+AU^  ,Vi+AV^  ) . 

For  this  example  definition  (5-22)  can  again  lead  to  erroneous  results. 
Again  the  assumed  principal  value  will  affect  the  average.  That  Is, 
whether  the  vector  In  quadrant  III  of  Figure  5.31c  Is  considered  to  be 
In  the  range  ir  to  3it/2  or  In  the  range  -it/2  to  -it  will  affect  the 
phase  average.  Also,  whether  the  vector  In  quadrant  IV  of  Figure  5.31c 
Is  considered  to  be  In  the  range  3ii/2  to  2it  or  In  the  range  0 to  -ii/2 
will  affect  the  phase  average. 

This  problem  may  be  circumvented  as  follows.  The  complex  plane  Is 
divided  Into  phase  zones.  For  example.  In  Figure  5.32  the  complex 
plane  has  been  divided  Into  eight  phase  zones.  A histogram  can  be 
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Figure  5.31  a)  Representation  of  phase  difference  corresponding  to 

OIF  b)  Representation  of  phase  differences  corresponding 
to  undegraded  subimages  c)  Representation  of  phase 
differences  corresponding  to  degraded  subimages 
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constructed  of  the  number  of  phases  included  in  the  average  which  fall 
into  each  of  the  phase  zones.  Using  the  histogram  information,  the 
general  orientation  of  the  phase  average  can  be  estimated.  Denoting 
this  phase  by  0(u,v),  all  phases  in  the  average  are  now  made  relative  to 
0(u,v).  The  relative  phases  are  averaged.  This  average  is  then  added 
to  e(u,v)  for  the  final  average. 

For  example,  with  the  eight  phase  zones  illustrated  in  Figure  5.32 
the  histogram  of  the  phases  illustrated  in  Figure  5.30  would  indicate 
that  the  general  orientation  of  the  phase  average  is  around  w.  Relative 
to  IT,  phases  counterclockwise  to  n would  be  considered  as  positive  and 
phases  clockwise  to  n would  be  considered  negative.  The  average  of  the 
relative  phases  would  be  around  0,  and  the  final  average  would  be  n + 0 = 
n. 

All  phase  averaging  within  this  research  used  the  above  technigue. 

5. 7 On  the  Relaxation  of  One  Assumption 

Because  of  the  difficulty  inherent  in  approximation  (5-8),  it  is 
perhaps  instructive  to  eliminate  this  approximation  . With  the  isolation 
and  elimination  of  the  problems  inherent  in  approximation  (5-8)  one  can 
study  the  system  performance  and  in  particular  the  adequacy  of  the  esti- 
mate of  the  phase  averages 

[0c  (u,v)  - 0p  (u+Au,v+Av)]  , 

•“i  M 

the  adequacy  of  the  phase  averaging  technique,  and  additionally  the 
adequacy  of  the  rates  of  convergence  of  the  phase  averaging. 

This  approximation  can  be  eliminated  in  the  simulations  by  the 

following  method.  Recall  that  within  a half  width  of  the  PSF  from  the 

116 


i edge  of  a subimage,  approximation  (5-8)  Is  not  approximate  but  equal. 

As  a result.  If  one  zeros  the  edges  of  the  subimages  of  the  undegraded 
I Image  to  an  extent  greater  than  a half  width  of  the  PSF  and  then  blurs 

i the  resulting  Image  via  convolution,  approximation  (5-8)  Is  now 

equality  everywhere. 

For  example.  Figure  5.33a  Is  an  example  of  an  Image  zeroed  at  the 
edges  of  the  subimages  and  Figure  5.33b  Is  the  result  after  blurring. 

The  estimate  of  the  average  phase  differences  of  the  undegraded 
Image  Is  no  longer  zero.  In  this  case  the  average  phase  differences 
will  converge  to  the  phase  differences  associated  with  the  subimage  of 
I Figure  5.34.  One  can  use  the  mathematical  function  such  as  Illustrated 

I In  Figure  5.34  to  estimate  the  average  phase  differences  associated 

« 

‘ with  the  undegraded  Image  or  one  can  use  a similar  zeroed  undegraded 

I Image  to  estimate  the  average  phase  differences  associated  with  the 

undegraded  image.  The  latter  approach  was  used  in  the  simulations  to 
' follow. 

If  one  zeros  the  edges  of  the  subimages  of  an  image  and  then 
calculates  the  average  phase  differences,  one  could  then  display  the 
I average  phase  differences.  On  the  other  hand,  a more  esthetically 

pleasing  display  is  usually  obtained  by  starting  with  the  phase  at  the 
DC  term  of  zero  and  using  the  average  phase  differences  to  recursively 
> calculate  a continuous  phase  plane. 

Accordingly,  zeroing  the  edges  of  the  subimages  corresponding  to 
Figure  5.1b,  calculating  the  associated  average  phase  differences, 
using  the  phase  of  the  DC  term  of  zero  and  recursively  calculating  the 
phase  plane  from  the  average  phase  differences,  the  phase  plane  dis- 
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Figure  5.33  An  image  zeroed  at  the  edges  of  the  subimages 
and  then  blurred 
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played  In  Figure  5.35a  was  obtained.  This  Is  the  same  phase  plane  as 
would  be  obtained  by  Fourier  transforming  the  function  In  Figure  5.34 
and  displaying  the  unwrapped  phase  plane.  Plotted  In  Figure  5.35b  Is 
the  phase  along  the  u-axis  of  the  phase  plane,  while  plotted  In 
Figure  5.35c  Is  the  phase  along  the  v-axis  of  the  phase  .plane. 

5 . 8 Suninary 

Recall  that  the  success  or  failure  of  the  recursive  estimates  of 
the  magnitude  and  phase  of  the  OTF  for  real  world  blurred  Imagery  Is 
contingent  on  at  least  five  conditions:  1)  whether  or  not  the  mathe- 
matical model  Is  adequate,  2)  whether  or  not  the  unknown  quantities 
In  the  estimates  can  be  adequately  estimated,  3)  whether  or  not  the 
averages  converge  adequately,  4)  whether  or  not  the  Iterations  are 
stable,  and  5)  whether  or  not  the  error  Inherent  In  the  approximations 
are  acceptable.  Additionally,  recall  the  success  or  failure  of  the 
phase  recursion  is  dependent  on  the  adequacy  of  techniques  which 
assume  only  knowledge  of  the  principal  value  of  the  phase. 

Table  5.4  summarizes  the  system  performance  with  respect  to  the 
above  contingencies. 

The  ill-conditioned  nature  of  the  phase  causes  most  approximations 
to  be  suspect.  In  particular,  the  approximation 

eg^(u,v)  0^(u,v)  + 0p^(u,v)  (5-23) 

is  not  a particularly  good  approximation  In  general.  Windowing  can 
sharpen  the  approximation,  but  at  the  same  time  the  windowing  causes 
a loss  of  resolution  in  the  estimate  of  the  phase  differences.  In  the 
simulations  presented  In  Chapter  6 it  will  be  demonstrated  that  with 
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•e  5.35  Phase  plane 
zeroed  at  t 


PHASE  ON  V-flXIS  PHASE  ON  U-AXIS 


Magnitude 

recursion 

Phase 

recursion 

Model 

adequate 

adequate 

Estimate  of  unknown 
quantities 

adequate 

adequate 

Convergence  of 
averages 

adequate 

adequate 

Stability  of 
iteration 

conditionally 

unstable 

stable 

Approximations 

acceptable 

acceptable 

conditionally 

unacceptable 

Adequacy  of 
assuming  only 
principal  value 
of  the  phase 

adequate 

Table  5.4  Sunmary  of  adequacy  of  contingencies  upon  which 
recursion  depends 
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the  elimination  of  this  approximation,  via  the  method  of  section  5.7, 
the  phase  recursion  gives  reasonable  estimates  of  the  phase. 

Windowing  is  one  approach  to  sharpening  the  approximation.  An- 
other approach  is  to  increase  the  subimage  size.  For  example,  if  the 
PSF  is  nonzero  over  an  extent  of  7 x 7 pixels, for  a subimage  size  of 
64  X 64  pixels  the  relationship 

g^(x,y)  ^ b(x,y)*f|x,y)  (5-24) 

will  be  exact  except  for  at  most  19%  of  the  64  x 64  pixels.  If  the 
subimage  size  is  256  x 256  pixels,  relationship  (5-24)  will  be  exact 
except  for  at  most  5%  of  the  256  x 256  pixels.  At  some  point  the  sub- 
image size  will  be  large  enough  such  that  the  approximation  (5-23)  will 
contain  sufficient  accuracy  to  support  a reasonable  estimate  of  the 
phase. 

Note,  however,  that  given  an  image,  as  the  subimage  size  in- 
creases the  number  of  subimages  of  the  given  image  decreases.  Thus,  it 
may  be  necessary  to  have  a very  large  image  and  very  large  subimages 
for  approximation  (5-23)  to  support  an  accurate  estimate  of  the  phase. 


Chapter  6 

RESULTS  OF  SIMULATIONS 
6.1  Presentation  of  Results 

In  this  chapter  will  be  presented  the  results  of  computer  simula- 
tions for  a variety  of  blurs.  The  format  of  presentation  will  be  as 
follows.  First,  a diagram  of  the  degrading  PSF  will  be  presented.  This 
will  be  followed  by  a comparison  of  the  magnitude  of  the  OTF  and  the 
corresponding  estimate.  Because  of  the  unstable  nature  of  the  recursive 
estimate  of  the  magnitude  of  the  OTF,  the  magnitude  estimate  was  cal- 
culated via  the  method  of  Cannon  [26-28].  Next,  comparisons  of  the 
phase  of  the  OTF  against  two  estimates  of  the  phase  will  be  presented. 

The  first  phase  estimate  was  obtained  using  recursion  (5-2)  and  the 
method  outlined  in  section  5.7;  the  second  phase  estimate  was  obtained 
using  Parzen  windowing  and  using  recursion  (5-13).  Recall  the  method 
of  section  5.7  relaxed  one  assumption.  The  first  estimate  will  be 
referred  to  as  Estimate  1 and  the  second  estimate  will  be  referred  to  as 
Estimate  2. 

The  degraded  image  together  with  seven  restorations  will  then  be 
presented.  The  assumptions  underlying  the  seven  restorations  are 
given  in  Table  6.1. 

As  regards  the  diagram  of  the  degrading  point  spread  functions, 

■ « 

the  point  spread  functions  are  assumed  to  be  lossless.  Thus,  the 
volumes  under  the  point  spread  functions  are  assumed  to  be  unity. 

All  plots  comparing  the  magnitude  of  the  OTF  are  of  the  log^Q 
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given 
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1 c 

given 

estimated  as  0 
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given 

Estimate  2 
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estimated 

given 
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estimated 

estimated  as  0 

estimated 

j!;  9 

estimated 

Estimate  2 

estimated 

estimated 

Estimate  1 

estimated 

Table  6.1  Key  to  restorations 
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of  the  magnitude  and  are  all  presented  on  the  same  scale.  For  compari- 
son purposes  values  of  the  log-ig  of  the  magnitude  be1ow-2  have  been 
set  to  -2  and  in  all  cases  the  bottom  of  the  skirt  is  at  the  -2.1 
level.  Because  the  point  spread  functions  are  assumed  to  be  loss- 
less, the  log^Q  of  the  DC  term  of  the  point  spread  function  will  be  0. 
Accordingly,  the  log^g  magnitude  plot  of  the  optical  transfer  functions 
will  range  from  0 to  -2.1,  the  bottom  of  the  skirt. 

The  estimate  of  the  power  spectrum  of  the  undegraded  image,  i.e. 
of  the  image  of  Figure  5.1a,  was  the  average  power  spectrum  of  the 
power  spectra  of  jhe  Images  of  Figures  5.1b,  c,  and  d. 

Phase  Estimate  1 used  the  image  of  Figure  5.1b  as  outlined  in 
section  5.7  to  estimate  the  phase  differences  associated  with  the 
undegraded  image.  In  addition,  phase  Estimate  1 used  subimages  of 
64  X 64  pixels  and  no  overlapping  of  subimages.  This  resulted  in  an 
averaging  over  64  subimages. 
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Phase  Estimate  2 estimated  the  average  phase  differences  corres- 
ponding to  the  undegraded  image  as  zero.  As  did  phase  Estimate  1, 
phase  Estimate  2 used  a subimage  size  of  64  x 64  pixels.  Unlike  phase 
Estimate  1,  however,  phase  Estimate  2 used  a S0%  overlapping  of  sub- 
images. Thus,  the  averaging  was  over  225  subimages. 

In  each  figure  which  compares  the  phase  of  the  OTF  and  the  corres- 
ponding estimates,  the  plot  scaling  is  identical  within  the  given  figure; 
on  the  other  hand,  to  adequately  present  different  ranges  of  phase 
variations  for  different  point  spread  functions,  the  scaling  is  usually 
different  from  figure  to  figure. 

6.2  Discussion 

Upon  comparison  of  the  magnitudes  of  the  optical  transfer  functions 
and  the  estimates  of  the  magnitudes,  it  is  clear  that  the  estimates 
adequately  retain  the  essential  features.  The  magnitude  estimates 
resolve  the  essential  features  to  at  most  1 1/2  orders  of  magnitude. 

Most  would  agree  that  this  is  fairly  good,  especially  in  light  of  the 
fact  that  the  degraded  images  which  were  input  to  the  system  were  given 
to  at  most -8  binary  digits  (bits)  of  accuracy. 

Upon  comparison  of  the  phases  of  the  optical  transfer  functions 
and  the  estimates  of  the  phases,  it  is  evident  that  phase  Estimate  1 
is  reasonable  near  the  axes  and  not  as  accurate  at  a distance  from 
the  axes.  Because  images  in  general  have  little  energy  concentration 
in  those  frequency  components  away  from  the  axes,  the  poorer  results 
away  from  the  axes  are  not  very  important  from  a restoration  point  of 
view. 
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Some  of  the  differences  observed  between  the  phases  of  the  optical 
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transfer  functions  and  phase  Estimates  1 are  due  to  the  unwrapping. 

It  should  be  emphasized  that  a consistent  unwrapping  is  necessary  for 
an  esthetically  pleasing  display,  but  not  for  the  Wiener  filter.  Some 
unwrapping  problems  are  evident.  For  example,  in  Figure  6.4a  it  may  be 
noted  that  at  one  point  the  display  of  the  phase  of  the  OTF  jumped  -n 
and  at  the  same  point  the  display  of  phase  Estimate  1 jumped  +tt.  Other 
unwrapping  problems  are  apparent,  for  example,  see  Figures  6.14b,  6.19a, 
and  6.19b. 

In  addition,  some  of  the  differences  observed  between  the  phases 
of  the  optical  transfer  functions  and  phase  Estimates  1 are  due  to 
the  tendency  of  the  phase  estimates  to  contain  a slight  linear  term. 

The  slight  linear  term  causes  only  a slight  shift  in  the  picture  domain 
and  is  not  responsible  for  any  blurring;  thus,  it  is  of  minor  conse- 
quence. However,  on  comparing  the  displays  of  the  phases  of  the  optical 
transfer  functions  and  the  displays  of  the  corresponding  estimates,  the 
comparisons  are  less  favorable.  For  example,  upon  comparing  the  phase 
of  the  OTF  and  Estimate  1 in  Figures  6.9a  and  6.9b,  it  is  evident  that, 
if  not  for  a slight  linear  term  in  Estimate  1,  the  estimate  would  be 
even  closer.  This  effect  is  also  evident  in  Figures  6.19a  and  b, 

6.24b,  and  6.29  a and  b. 

From  the  phase  comparisons  it  may  be  noted  that  phase  Estimate  2 
is  not  particularly  good. 

Comparing  restorations,  it  is  apparent  that  the  estimate  of  the 
magnitude  of  the  OTF  is  reasonable  to  the  extent  of  achieving  an 
adequate  restoration.  In  fact,  the  visual  quality  of  the  restorations 

using  the  estimate  of  the  magnitude  of  the  OTF  seems  to  be  superior  to 
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the  restorations  given  knowledge  of  the  magnitude  of  the  OTF.  For 
example,  compare  Figures  6.10  b and  e.  Figures  6.15  b and  e.  Figures 
6.20  b and  e,  and  Figures  6.30  b and  e.  In  general,  upon  restoration, 
the  frequencies  in  the  vicinity  of  the  zeros  of  the  OTF  will  be 
magnified  the  most  by  the  Wiener  filter.  If  one  knows  the  magnitude 
of  the  OTF  exactly,  at  times  the  restoration  will  contain  periodic 
artifacts  corresponding  to  those  frequencies  that  were  boosted  the  most 
by  the  filter.  On  the  other  hand,  the  inaccuracies  inherent  in  the 
estimate  of  the  magnitude  of  the  OTF  cause  the  frequencies  in  the 
vicinity  of  the  zeros  to  be  boosted  less  than  the  boost  assuming 
knowledge  of  the  magnitude  of  the  OTF.  As  a result,  when  assuming 
the  estimate  of  the  magnitude  of  the  OTF,  in  most  cases  the  periodic 
artifacts  are  missing. 

Although  the  visual  quality  of  the  restoration  is  usually  superior 
when  using  the  estimate  of  the  magnitude  of  the  OTF,  in  terms  of 
minimum  mean  squared  error,  restoration  b is  superior  in  all  cases. 

Upon  comparing  restorations  c and  d of  a set  and  restorations  f 
and  g of  a set,  it  is  generally  agreed  that  Estimate  2 of  the  phase  is 
not  significantly  better  with  respect  to  restoring  than  an  estimate  of 
zero  for  the  phase. 

The  superiority  of  Estimate  1 of  the  phase  is  evident.  For  example, 
note  in  Figures  6.5  f and  g the  double  set  of  eyes  and  the  collar  of  the 
man;  also  note  the  streak  of  light  across  the  arm  and  chest  of  the  t 

woman  and  below  the  knees  of  the  womaa  These  anomalies  are  absent  in 
Figure  6.5h.  Note  in  Figure  6.15h  the  elimination  of  the  double  image 
in  Figures  6.15  f and  g.  Note  the  clearer  definition  of  the  faces  in 
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Figure  6.25h  opposed  to  the  faces  In  Figures  6.25  f and  g.  Additionally, 
note  the  general  superiority  of  Figure  6.30h  over  Figures  6.30  f and  g. 
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(a)  degraded  image 


magnitude  of  OTF  given 
Phase  estimate  = Estimate  2 


Figure  6.5  Restorations 


(e)  magnitude  of  OTF  estimated  (f)  magnitude  of  OTF  esti- 

phase  of  OTF  given  mated,  phase  of  OTF  esti- 

mated as  0 


(g)  magnitude  of  OTF  estimated  (h)  magnitude  of  OTF  esti- 

phase  estimate  = Estimate  2 mated, phase  estimate  = 

Estimate  1 
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Figure  6.5  (Cont'd) 
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(b)  V axis  results 

Figure  6.9  Comparison  of  phase  of  OTF  and  estimates 


(a)  degraded  image 


(b)  magnitude  of  OTF  given 
phase  of  OTF  given 


(c)  magnitude  of  OTF  given 

phase  of  OTF  estimated  as  0 


(d)  magnitude  of  OTF  given 
phase  estimate  = 
Estimate  2 


Figure  6.10  Restorations 
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(e)  magnitude  of  OTF  estimated 
phase  of  OTF  given 


(f)  magnitude  of  OTF  esti- 
mated, phase  of  OTF 
estimated  as  0 


(h)  magnitude  of  OTF  esti- 
mated, phase  estimate  = 


^g)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  2 
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Fiuqre  6.10  (Cont'd) 
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(a)  u axis  results 
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(b)  V axis  results 

Figure  6.14  Comparison  of  phase  of  OTF  and  estimates 


(a)  degraded  image 


(b)  magnitude  of  OTF  given 
phase  of  OTF  given 


(d)  magnitude  of  OTF  given 

phase  estimate  = Estimate  2 


Figure  6.15  Restorations 
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magnitude  of  OTF  given 
phase  of  OTF  estimated  as  0 
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(e)  niaqnitude  of  OTF  estimated 
phase  of  OTF  given 


(f)  magnitude  of  OTF  estimated 
pha^e  of  OTF  estimated  as  0 


(q)  magnitude  of  OTF  estimated  (h)  magnitude  of  OTF  estimated 

phase  estimate  = Estimate  2 phase  estimate  - Estimate  1 


Figure  6.15  (Cont'd) 
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Estimate  1 


Phase  of  OTF 


Estimate  2 


FREQUENCY 


u axis  results 


Phase  of  OTF 


stimate  1 


Estimate  2 


FREQUENCY  ^ 

(b)  V axis  results 

Figure  6.19  Comparison  of  phase  of  OTF  and  estimates 
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(a)  degraded  image 


(b)  magnitude  of  OTF  given 
phase  of  OTF  given 


(c)  magnitude  of  OTF  given 

phase  of  OTF  estimated  as  0 


(d)  magnitude  of  OTF  given 

phase  estimate  = Estimate  2 


Figure  6.20  Restorations 
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(e)  magnitude  of  OTF  estimated 
phase  of  OTF  given 


(f)  magnitude  of  OTF  estimated 
phase  of  OTF  estimated  as  0 


(g)  mag.-itude  of  OTF  estimated 
phase  estimate  = Estimate  2 


(h)  magnitude  of  OlF  estimated 
phase  estimate  = Estimate  1 


Figure  6.20  (Cont'd) 
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Estimate  1 


Phase  of  OTP 


Estimate  2< 


Figure  6.2 
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(e)  magnitude  of  OTF  estimated 
phase  of  OTF  given 


(f)  magnitude  of  OTF  estimated 
phase  of  OTF  estimated  as  0 


(g)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  2 


(h)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  1 


Figure  6.25  (Cont'd) 
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(a)  e 


(b)  0=Estimate  1 


(c)  0*Estimate  2 

Figure  6.28  Comparison  of  phase  of  OTF  and  estimates 
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(a)  degraded  image  (b)  magnitude  of  OTF  given 


phase  of  OTF  given 


(c)  magnitude  of  OTF  given  (d)  magnitude  of  OTF  given 

phase  of  OTF  estimated  as  0 phase  estimate  = Estimate  2 


Figure  6.30  Restorations 
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(e)  magnitude  of  OTF  estimated 
phase  of  OTF  given 


(f)  magnitude  of  OTF  estimated 
phase  of  OTF  estimated  as  0 


(g)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  2 


(h)  magnitude  of  OTF  estimated 
phase  estimate  = Estimate  1 


Figure  6 30  (Cont’d) 
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Chapter  7 

RESULTS  ON  REAL  WORLD  BLURRED  IMAGES 

This  chapter  presents  the  results  for  three  photographically 
Induced  blurs.  To  Induce  the  first  two  blurs  the  camera  was  jiggled 
and  vibrated  during  exposure.  The  third  blurred  image  was  obtained 
from  a private  source;  the  blur  was  apparently  the  result  of  an 
incorrect  use  of  the  camera. 

The  blurred  photographs  were  digitized  to  512  x 512  pixels  and 
the  estimates  of  the  OTF  were  made  using  the  512  x 512  blurred  images. 

A 50X  overlapping  of  subimages  gave  225  subimages  of  64  x 64  pixels 
each.  The  estimate  of  the  magnitude  of  the  OTF  was  made  via  the 
method  of  Cannon  [26-28].  The  estimate  of  the  power  spectrum  of  the 
undegraded  image  was  the  average  power  spectrum  of  the  power  spectra 
of  the  images  of  Figures  5.1b,  c,and  d.  The  phase  estimate  was  obtained 
using  a Parzen  window  and  using  recursion  (5-13). 

For  the  first  two  blurred  images  the  restorations  were  made  on 
each  of  the  four  256  x 256  pixel  quadrants  of  the  512  x 512  pixel 
blurred  image.  For  the  third  blurred  image  four  restorations  of 
256  X 256  pixels  were  made;  the  four  restorations  centered  prominent 
features  within  the  512  x 512  pixel  blurred  image. 

In  Figure  7.1  is  the  same  scene  before  and  after  the  photographic- 
ally induced  blur.  Note  the  "before"  picture  is  earlier  in  time  and 
is  not  used  in  the  restoration.  Figure  7.2  presents  perspective  plots 
of  the  estimates  of  the  magnitude  and  phase  of  the  OTF. 

Figures  7. 3-7. 6 present  the  results  of  the  restorations  on  the 
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Figure  7.1  Scene  before  and  after  photographically  induced  blur 
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(b)  e (minimum  * -1.04  radians, 
maximum  > 1.04  radians) 

Figure  7.2  Estimates  of  magnitude  and  phase  of  OTF 


(a)  degraded 


(b)  phase  of  OTF  estimated  as  0 (c)  phase  estimate  = Estimate  2 


Figure  7.3  Blurred  image  and  restorations 


1 70 


(b)  phase  of  OTF  estimated  as  0 (c)  phase  estimate  = Fstimate  ? 


Figure  7.4  Blurred  image  and  restorations 
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(a)  degraded 


(b)  phase  of  OTF  estimated  as  0 (c)  phase  estimate  = Estimate  2 


Figure  7.5  Blurred  image  and  restorations 
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(a)  degraded 


(b)  phase  of  OTF  estimated  as  0 (c)  phase  estimate  = Fstimate  2 


Figure  7.6  Blurred  image  and  restorations 
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four  256  X 256  pixel  quadrants  of  the  512  x 512  pixel  blurred  image. 

In  each  of  Figures  7. 3-7.6  is  presented  the  degraded  quadrant  together 
with  a restoration  using  the  estimate  of  the  magnitude  and  phase  of 
the  OTF  and  a restoration  using  the  estimate  of  the  magnitude  of  the 
OTF  and  estimating  the  phase  of  the  OTF  to  be  zero. 

Improvement  is  evident.  In  addition  to  an  improvement  in  sharp- 
ness, some  objects  that  were  not  recognizable  in  the  blurred  images 
are  recognizable  after  restoration.  In  Figure  7.3a  the  vertical 
columns  of  the  building  in  the  northeast  quadrant  of  the  image  are  not 
resolved.  In  the  restorations  in  Figures  7.3b  and  c the  columns  are 
resolved. 

In  Figure  7.4a  the  object  on  top  of  the  tower  is  not  recognizable. 
In  the  restorations  in  Figures  7.4b  and  c it  is  seen  to  be  a ball.  In 
addition,  in  Figures  7.4b  and  c note  the  improvement  in  definition  of 
the  windows  and  structure  of  the  building  which  is  in  the  center  of 
the  right-hand  side  of  the  frame. 

In  the  restorations  of  Figures  7.5b  and  c the  lines  of  the 
crosswalk  are  now  defined.  Additionally,  there  is  better  definition  in 
the  cars;  it  is  now  possible  to  recognize  the  Volkswagen  as  a Volkswagen. 
Note  the  increased  resolution  in  the  windows  of  the  minibus. 

In  Figures  7.7-7.12  are  presented  the  results  corresponding  to 
the  second  photographically  induced  blur.  Again,  an  improvement  in 
sharpness  arJ  increased  resolution  is  observed.  In  Figures  7.9b  and  c 
note  the  increased  sharpness  and  definition  in  the  tree  in  the  south- 
west quadrants,  the  tree  in  the  center  of  the  right-hand  side  of  the 

frames,  and  the  trees  along  the  top  of  the  frames.  Additionally,  note 
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(b) 

Figure  7.7  Scene  before  and  after  photographically  induced  blur 
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(b)  phase  of  OTF  estimated  as  0 (c)  phase  estimate  = Estimate  2 

Figure  7.9  Blurred  image  and  restorations 
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(b'  phase  of  OTF  estimated  as  0 


(c)  phase  estimate  = Estimate  2 


Figure  7.10 


Blurred  image  and  restorations 
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(a)  degraded 


(b) 


phase  of  OTF  estimated  as  0 (c)  phase  estimate 


Figure  '.11  Blurred  image  and  restorations 


Estimate  2 
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(a)  degraded 


(b)  phase  of  OTF  estimated  as  0 


(c)  phase  estimate  = Estimate  2 


Figure  7.12  Blurred  image  and  restorations 
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the  increased  sharpness  in  the  cars  and  building. 

In  Figures  7.10b  and  c there  is  better  definition  in  the  window 
panes  and  ledges.  Note  the  increased  sharpness  in  the  trees  of 
Figures  7.12b  and  c. 

Figures  7.13-7.15  present  the  results  for  the  third  real  world 
blurred  image.  Figure  7.13  presents  the  estimate  of  the  magnitude  and 
phase  of  the  OTF.  Figures  7.14  and  7.15  present  restorations  of 
selected  256  x 256  pixel  subimages  of  the  larger  512  x 512  blurred 
image.  Again,  improvement  is  evident.  For  example,  in  Figure  7.14b 
note  the  increased  sharpness  of  the  edge  content  compared  to  Figure 
7.14a.  Additionally,  in  Figure  7.14b  note  the  better  definition  in 
the  rocks. 

In  Figures  7.14d  and  7.15b  in  addition  to  an  increased  sharpness 
some  letters  that  were  not  discernible  in  the  blurrev^  versions  are  now 
discernible. 

In  Figure  7.1b  there  are  what  appears  to  be  degraded  point 
sources.  For  example,  note  what  appears  to  be  a degraded  point  source 
slightly  below  the  center  of  the  frame  on  the  left-hand  side  of  the 
frame.  This  affords  one  with  the  opportunity  of  estimating  the 
magnitude  and  phase  of  the  OTF  using  the  degraded  point  source. 

In  Figure  7.16  is  a comparison  of  the  results  using  the  techniques 
forming  the  basis  of  the  research  reported  herein  and  results  obtained 
using  a degraded  point  source.  Figures  7.16a  and  b are  the  same  as 
Figures  7.5a  and  c,  respectively,  repeated  here  for  convenience. 

Figures  7.16a  is  the  degraded  Image.  Figure  7.16b  is  the  restoration 
obtained  using  the  techniques  reported  herein. 
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Figure  7.13  Estimates  of  magnitude  and  phase  of  OTF  for  third 
real  world  blurred  Image 
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(a)  degraded 


(c)  degraded 


Figure  7.14  Blurred  images 


(d)  restored 


snd  restorations 
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(a)  degraded 


(b)  restored 


(c)  degraded 


(d)  restored 


Figure  7.15  Blurred  images  and  restorations 
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(a)  degraded 


(b)  restoration  using 

techniques  reported  herein 


(c)  restoration  using  degraded 
point  source 


(d)  magnitude  of  OTF  via  method 
of  Cannon,  phase  of  OTF 
from  degraded  point  source 


Figure  7.16  Comparison  of  restoration  results  with 
restoration  results  using  the  degraded 
point  source 
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Using  the  degraded  point  source  In  the  northwest  quadrant  of 
Figure  7.16a  to  estimate  both  magnitude  and  phase  of  the  OTF  and 
using  these  estimates  to  estimate  the  undegraded  Image,  Figure  7.16c 
was  obtained.  The  degraded  point  sourcesare  now  closer  to  points. 
However,  the  rest  of  the  Image  Is  probably  a little  worse.  It  Is 
concluded  that,  although  theoretically  possible.  It  Is  very  difficult 
In  practice  to  restore  blurred  Images  by  using  degraded  point  sources 
to  estimate  the  OTF. 

Figure  7.16d  uses  the  magnitude  estimate  used  In  Figure  7.16b, 
that  Is,  via  the  method  of  Cannon  [26-28],  and  uses  the  phase  estimate 
from  the  degraded  point  source.  Note  the  degraded  point  sources  are 
closer  to  points.  Note,  additionally,  the  improvement  in  the  definition 
of  the  back  of  the  Pinto  parked  Immediately  before  the  crosswalk. 
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Chapter  8 


SUMMARY 


8.1  Summar 


Because  of  the  unstable  nature  of  the  recursive  estimate  of  the 


magnitude  of  the  OTF,  this  estimate  is  of  little  practical  value 


With  regards  to  the  recursive  estimate  of  the  phase  of  the  OTF, 
from  the  results  of  Chapter  6 it  it  seen  that,  via  the  method  of 
section  5.7,  a reasonable  estimate  of  the  phase  of  the  OTF  can  be 
calculated.  Recall  the  method  of  section  5.7  zeros  the  edges  of  the 
subimages  prior  to  the  blurring  operation.  This  has  the  associated 
effect  of  sharpening  the  approximation 


Another  approach  to  sharpening  the  approximation  is  to  use  larger 
subimages  in  the  calculations.  As  the  subimage  size  is  increased, 
however,  the  cost  of  the  phase  estimate  increases  accordingly.  In 
addition,  a larger  image  is  needed  to  support  a reasonable  convergence 
of  the  averages. 


8.2  Open  Question 


Given  a PSF  of  a given  nonzero  extent,  at  what  subimage  size  will 


the  approximation 


APPENDIX  A 

Theorem:  The  trigonometric  polynomial 

T(u)  . 2 
k-O 

can  have  at  most  2L  real  2eros  In  the  Interval  j j . even  If  each 
multiple  root  Is  counted  the  number  of  times  It  occurs. 

Proof:  Let  us  use  the  following  Identities. 

cos  u 

It  follows  that 

T(u) 

where  the  C|^'s  are  complex  constants. 

T(u)  • ^ c^j.(j2i«(k+L)u/N) 

k— L 

. g(-J2nLu/N)  ^ j^g(j2itiu/N) 
l-O 

Now  consider  the  polynomial 
P(Z)  ■ 


gju  ^ g-Ju  e''”  - e 


Ju  _ -ju 
, sin  u • jj- 


Z 'k' 

k— L 


(j2nku/N) 


and  form  the  correspondence 
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j . e(J2i.u/N) 

Note  that  to  each  uc  ^ j corresponds  a unique  z.  Now 

p(z)  - . 

Differentiation  with  respect  to  u gives 

3P(z)  ii  . iliiL  e(J2wLu/N),,  » . ^(,i2nLu/N)  8T(u) 
32  au  N ® * ® ~3u"^ 


and  via  Induction  It  can  be  shown  that 


p(s)(2)  . [XqT(u)  + a^T'(u)+...+  Xj.T^*^(u)]  (A-1) 

where  the  x^'s  are  complex  constants. 

Now  assume  that  Uq  Is  an  m-fold  zero  of  T(u).  This  Implies  that 

T(uq)  - T^Uq)  T^'"“^^(uq)  » 0. 

For  Zq  ■ equation  (A-1)  It  follows  that 

1 

P(Zq)  ■ p^Zq)  p^^’^^z^)  - 0. 

Thus,  to  each  m-fold  zero  of  T(u)  corresponds  an  m-fold  zero  of  p(z). 

The  nunberof  zeros  of  p(z)  cannot  exceed  2L  where  the  multiplicities 
of  each  zero  Is  taken  Into  account.  Since  to  each  zero  of  T(u)  of 
multiplicity  m for  uc  f-  ^ ,^j  corresponds  d zero  of  p(z)  of 
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